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ABSTRACT 


The waves propagating in media with transition 
layers are investigated theoretically with application to 
the propagation of seismic waves in the Earth's crust and 
the upper mantle. The propagation of pressure perturbation 
from a point source in a fluid is studied in both the frequency 
and time domain. The problem is solved in the flat geometry 
under the assumption that the earth-flattening transformation 
can be used to transform the spherical geometry into the 
flat geometry. 

When the monotonic velocity transition is broad the 
caustics are formed by the focussing of rays due to the rapid 
change in the velocity gradient. In such a case a triplication 
is formed on the travel-time curve and the wave propagation in 
the neighbourhood of the endpoints is studied. Numerical 
evaluation of the response of the medium shows that at high 
frequencies the decay into the region beyond the caustic is 
large and that low frequency waves carry the significant 
amount of energy. The amplitude maximum is shifted from the 
endpoint into the eho Ant cuca region due to constructive 
interference of the two arriving waves. 

If the monotonic transition is thin the partial 
reflection occurs at sub-critical angles of incidence and the 


head wave arises at the critical incidence. They are both 
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effected by the thickness of the transition: as the transition 
becomes thinner, they behave more like sub-critical reflection 
and head waves from sharp velocity transition. 

When the transition layer has a velocity maximum 
the amplitudes of the waves propagating along the velocity 
reversal decay exponentially with distance. This rapid 
decrease of the amplitude causes a "shadow" whose properties 
are very different from a "true" shadow due to an abrupt 
velocity decrease. The first is frequency independent while 
the other depends on frequency as exp (aan) : 

The numerical results for the monotonic velocity 
transition are compared to the results given by the complete 
and partial ray expansions for medium with homogeneous layers. 
The results suggest that if a good approximation of the 
inhomogeneous medium by a stack of homogeneous layers is 
required the number of layers must be very large. The error 
of the partial ray expansion with respect to the exact solution 
depends on the choice of rays considered. For large number 
of layers the rays with multiple reflections must be included 


in the expansion in order to obtain satisfactory results. 
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CHAPTER 1 INTRODUCTION 


Improved coverage by stations of the USCGS world- 
wide network has provided better data in the last decade. 
Their processing has lead to more accurate estimates of elastic 
wave velocities and density within the Earth. All the recent 
models (Johnson, 1967; Julian and Anderson, 1968, for example) 
show remarkable agreement with the classical models of 
Jeffreys and Bullen (1940) and Gutenberg (1958) except at 
regions of high velocity gradients or low velocity channels 
(Figure 1.1). All these anomalous regions in the velocity 
structure are characterised by discontinuities or reversals 
in the travel-time curves. Increased resolution of these 
curves is the basis for better interpretation and leads to 
better understanding of the Earth's interior. This thesis 
studies theoretically a number of phenomena displayed by 
waves propagating in media with velocity transitions. 

The first models of the velocity structure within 
the Earth were derived by simple geometrical ray theory 
(Bullen, 1963; p. 109). It solves the problem of wave propagation 
when energy propagates within a ray tube without leaking 
through its walls. The vector of energy flux is normal to 
the wave front. The energy propagation cannot be explained 
so simply for many signals. For example the high velocity 


gradients in the upper mantle give rise to triplications on 
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Velocity structure of 
compressional waves in the 
upper mantle (from Johnson, 


1967.).. 
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the travel-time curves at approximately 15° and 20°. The 
later arrivals were. first carefully studied by Johnson (1967) 
and their interpretation lead to model CIT 204 Pier loel¢..tne 
ray tube on the caustic has a zero cross-section caused by 
focussing of the rays. Therefore the simple assumption about 
propagation of energy fails. The wave propagation in the 
region of the caustic was first studied by Jeffreys (1939) 
in connection with the PKP caustic. Other signals for which 
the geometrical ray theory is inapplicable are head waves in 
the vicinity of the critical angle. These were first studied 
by Cagniard (1939) and recently Cerveny and Ravindra (1971) 
published an excellent book that studies head waves in great 
detail. Diffraction signals cannot be explained by geometri- 
cal ray theory either (Knopoff and Gilbert, 1961; Alexander 
and Phinney, 1966; Phinney and Cathles, 1969 and others). 
In all these cases the wave equations that describe the 
motion of the medium have to be solved exactly in the frequency 
range needed as the response is frequency dependent. The 
geometrical ray theory forms the high frequency approximation. 
The solution must satisfy the source condition, the 
boundary conditions and the radiation conditions. This is 
equivalent to finding a reflection coefficient. The original 
equation is usually transformed to one variable. The trans- 


formation is in effect a decomposition of a spherical or 
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cylindrical wave into plane waves which are solutions of 
simpler differential equations. The wave equation for one 
variable must be solved and the reflection coefficient found. 
But even the one-dimensional equation can be solved analyti- 
cally for only a few models. For more complicated velocity 
structures numerical methods of solution must be employed 
(Chapman, 1969; Muller, 1970; Chapman and Phinney, 1972). 

With the great development of powerful computers the impor- 
tance of analytic solutions seems to have faded in recent 
years. The numerical solutions, however, carry with them a 
heavy computing penalty if we want to achieve the necessary 
accuracy. It seems useful to reinvestigate existing solutions 
and look to other fields of theoretical physics, if problems 
solved there are relevant to elastic wave propagation. Analytic 
solutions are still an indispensible part of research in 
theoretical seismology. However limited the model may appear 
with respect to the Earth's structure the results obtained 
often prove themselves to be very valuable. 

When the reflection coefficient has been found, the 
seismic response of any medium is obtained as a product of 
elementary waves that represent the source and the reflection 
coefficient. The complete solution which covers all frequen- 
cies can sometimes be obtained exactly in analytic closed form. 


That is the case when the reflection coefficient is independent 
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of frequency and the Cagniard-de Hoop method can be used 
(Cagniard, 1939; de Hoop, 1960). Gilbert and Helmberger 
(1972) expanded the spherical reflection coefficient into 

an asymptotic, frequency independent form in order to 

obtain solution for response of a layered sphere in time 
domain. Muller (1970) and Helmberger and Wiggins (1971) 
applied an earth-flattening approximation to the spherical 
model and then applied the exact Cagniard-de Hoop method to 
the approximate model of plane layers. But generally the 
reflection coefficient is frequency dependent and the response 
of the medium in the frequency domain must be found first. 

The time domain response can be obtained by inverse Fourier 
transformation. In many cases the results in frequency domain 
are more valuable as seismological data are usually analysed 
in the frequency domain. 

The small family of models for which the reflection 
coefficient has been found includes the homogeneous layers, 
linear transition layers (Gupta, 1966; Nakamura, 1964; 

Ward, 1972), exponential and parabolic bateok ted dee 197135 
Rydbeck, 1943) and Epstein layers (Epstein, 1930; Rawer, 1939; 
Lang and Shmoys, 1968; Phinney, 1970). Budden (1961) and 
Brekhovskikh (1960) summarize some of these results. Linear 
layers have discontinuities in velocity gradient which may 


not be realistic. Merzer's study is concerned only with head 
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waves and studies the variation of their amplitudes with 
frequency, layer thickness and shape of the transition. The 
Epstein models offer a wider range of applicability. We can 
model a region of high velocity gradient that gives rise to 
the triplication on the travel-time curve (Figure 1.2) or a 
narrow transition for investigation of the head wave (Figure 
1.3). We have already mentioned that the classical ray theory 
cannot describe the energy propagation in the neighbourhood 
of the caustics and around the critical point. Furthermore 
if the transition is continuous it cannot explain any sub- 
critical or partial reflection. Epstein (1930) and Phinney 
(1970) studied the modulus of the reflection coefficient 
to obtain an estimate of the amplitude of the partial reflec- 
tion. The modulus of the reflection coefficient itself, 
however, cannot give a good picture about energy propagation 
when the incident wave interacts with the strong velocity 
gradient. That can only be studied using the spectral 
amplitudes and impulse response of the model. A principle 
contribution of this work is their evaluation mae investiga- 
tion (Chapter 4). 

Another velocity variation we can study is the 
velocity reversal which does not cause any discontinuity on 
the travel-time curve. Nevertheless, at large ranges the 


signals carry very little energy due to large geometrical 
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Figure 1.2 The. triplication.and caustics due to the 
broad velocity transition. 
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Figure 1.3 The partial reflection and head wave from a 
marrow velocity transition. 
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spreading along the wave path. The low frequency signals leak 
through the velocity barrier. The "shadow" caused by velocity 
reversal has not been studied and the results can be useful 
for investigation of the upper part of the low velocity 
channel and the region D" at the base of the mantle. The 


"shadow" caused by this velocity structure is 


nature of 
basically very different from a true shadow formed by an 
interface with velocity decrease. The interface shadow has 
a frequency dependent diffracted signal (Scholte, 1956; 
Duwalo and Jacobs, 1959; Chapman, 1969; Chapman and Phinney, 


1972, and others) whereas the reversal ' 


"shadow" is indepen- 
dent of frequency. The discussion of the results for the 
reversal and the interface shadows is presented in Chapter 5. 
We have already mentioned that solution for wave 
propagation in inhomogeneous media can be obtained by numerical 
solution of the wave equation and that this requires sophisti- 
cated high-speed computers. Approximation of such media by 
a stack of homogeneous layers has often been used to solve the 
problem approximately. Solution for propagation of plane 
waves in a layered medium was found by Thomson (1950) and 
Haskell (1953). Cisternas et al. (1973) showed that such 
solution can be written as an infinite series where each 


term represents one ray, i.e. we can say that the exact 


solution is equivalent to the "complete" ray expansion. In 
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practise only a "partial" ray expansion is used, i.e. only 

the rays with greatest amplitude are considered in evaluation 
of the response of the medium (Muller, 1970; Hron, F., 1971). 
The rest of the series is neglected under the assumption that 
the series is convergent. The convergence has not been 

proved to date, and we will investigate it in Chapter 6. The 
estimate of the error introduced by the approximation of the 
inhomogeneous medium by stratified medium of homogeneous layers 
can be obtained by comparison of the complete ray expansion 
with the exact solution (analytic, when exists, or numerical). 
In spite of growing importance of partial ray expansion for 
theoretical seismology (Gilbert and Helmberger, 1972; Waddington, 
1973) very little is known about the error introduced by 

using it in place of the complete ray expansion. From this 
point of view the study presented in Chapter 6 can be regarded 


as a contribution to the solution of this complex problem. 
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CHAPTER 2 BASIC EQUATIONS 


ge ob Equations of motion 


The displacement equation of elastic motion may be 


found from the equation of conservation of momentum: 
sya, sito +p (x). £ (2.1) 


where u(x,t) is infinitesimal elastic displacement, 0(x) 
is the density and o (x) is the stress tensor. f£(x)is the 
applied force per unit mass acting at xX. Richards (1971) 


has introduced a potential representation of the displacement 


(r, 6,0) 


Ie 
ie) 
ct 

|x 

I 


which provides the correct separation in spherically inhomo- 
geneous media. The density o(r) and shear modulus u(r) are 
functions of radius and so is the scalar Function ¢ (rt); The 
potentials P(x,t), S(x,t) and T(x,t) can be considered as 
representing P,SV and SH wave motion in the high frequency 
limit. When (2.2) is substituted into Ce bd Be CE ee fg 
represented similarly, a second order differential equation 


for T(x,t) plus two coupled second order differential equations 
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for P(x,t) and S(x,t).are obtained. At high frequencies or 
when the medium is nearly homogeneous the coupling terms are 
small. Thus, at high frequencies we have three Helmholtz 

wave equations for P, SV and SH wave motion separately. Their 
solution can be used as approximate solutions for P and SV 
Waves and are exact for SH waves. We know from observations 
that the frequency range where P and SV waves are clearly 
distinguishable is quite wide. Thus the solutions of the 
Helmholtz equation are of great importance. 

We shall study motion in a fluid whose parameters 
vary smoothly with depth. We shall solve the problem by 
Epstein's method (Epstein, 1930) for pressure P rather than 
for potential Pp? - The density variation will be included 
in the solution. In Richard's formulation, the Helmholtz 


=e 
equation for potential Pp ” 


is valid only under the assumption 
of small density variations with respect to the wavenumber 
(Hill, 1971). This condition is violated for a discontinuity 
in density and, therefore, other methods of solutions must be 
employed in this case (Brekhovskikh, 1960; p. 171). The 
solutions for pressure P has no such disadvantage and is 
applicable to the decoupled P-SV motion or SH motion subject 


to the aformentioned constraints. 


Pressure in the fluid is given by 


Pit oel SH ee) ™ lx, t) (2.3) 
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and the equation of motion by 


Py aoe oe 8 Gh) (2.4) 


K(z) and 9(z) are coefficient of incompressibility and 
density, respectively. u(x,t) is the displacement at 
m(r,¢,2) and P(x,t) is the pressure. We assume that the 
source and the solution are azimuthally symmetric. In 


cylindrical coordinates (',¢,2) this means that i = 0 


ag : 
In practice this is not a serious limitation as the frequency 


range of the acoustic waves of interest to us is such that 


lateral inhomogeneities would effect the propagation very 
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if weMeliminate *a*="u(riopz,t) in.(2.4)%webget 
1 9p 1 
Sis geine beie Bt2ayv .* or Gz ¥ =» Vrp) (2.5) 
Vv ot 
Kz)... P : 
where v(z) = ate} is the velocity of acoustic pressure 


perturbations. Equation (2.5) is a second order differential 
equation and in order to solve such an equation analytically 


we must transform the variables. 


2.2 The coordinate transformatio 





First, equation (2.5) is Fourier transformed with 


respect to time, i.e. 
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+00 +00 +00 
f(t) = = if eon Wece vero dt' dw = a [ Fqwye7tet dw 
—©O —CO —co 


The Fourier transformation ig denoted by curly tilde. ‘f 
denotes any function of t. Further, the finite Fourier 


transform with respect to the angle Ob ia-uped.: 1.¢€ « 


27 
+0 : : 1 +o “a . 
£() -= gir 126 | fires. ad =-5 I £(2yei*? 
QR=—0 Q=-—00 
1@) 


ee a 


The finite Fourier transform is denoted by capped tilde. Since 
£(¢) must be a single valued function of >? , & must be an 
integer. In our case of azimuthal Symmetry £ = 0 is the only 
non-zero term. The Hankel transform or the Fourier-Bessel 


transform with respect to r is used, i.e. 


f(r) = if K Jo (kr) j vr" fieer' ) Jy (kr") dim’ tdic w= J K £(K) Jy (kr) dk 
re) re) ° 


(2.8) 


The Hankel transform is denoted by double capped tilde. 
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We also define the Laplace transformation of time 





+i 00 +ji0 
on on st j rw Bt) ; ed. st 
a a = ITd J e Ga ye dt ds = OTL i f(s) e ds 


(2.53 


When Laplace transformation is used to transform time t, 


the Laplace-Bessel transform is used to transform r. 


+ie 
2 
f(r) = 7 Im J P K, (spr) j al I,(spr') P(e} de* dp = 
-j{o“ fe) 
2 +j{° nm 
oe — Im j pr itp) Ko (spr) dp (2,403 
-j{o 


K, (spr) and I,(spr') are the modified Bessel functions. The 
equivalence of the Laplace-Bessel and Fourier-Bessel trans- 
formations is evident if we change the variables s = -iw and 
isp = K and convert the modified Bessel functions to I, (spr) 
and K) (spr) to Jo (kr). The relationship between the trans- 


formed functions is 


Dh per ftey, ss (2.11) 
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2.3 The sotirce 


Let us consider a point pressure source Pi) 
residing at r = 0, 2 = ze where the medium varies slowly. 


The homogeneous equation (2.5) is replaced by 


o(z) Ve (—L~ vp) 


p(z) St Te = Beers Ott) 6 Cane 2) (77103 


9*p 
str 


ig 
2 
Vv 


After transforming the equation (2.12) with respect to 


t, >, r we get the following equation for transformed 





pressure Pei= Bis 6 se) 

d 1 dP a 2 - 
P(z) a Cate ome a a - K )P = Pw) 6(z-z_) (2), 73) 
for £ = 0, and equal to zero for 2 #0. k= < is the wave- 
number, K = kesin 6 is the horizontal wavenumber’ and 

2 * 

Q = S - KK?) = kecos 6 is the vertical wavenumber. The 

Vv 


! 


Riemann surface Im Q > 0 is the "physical" Riemann sheet. 


This choice is bound to our choice of the Fourier transform. 
: : -iwt 
It assumes the time dependence of the wave being e 


1Wt ri : > 
rather than e which has often been used in earlier works 


on elastic wave propagation. Our choice makes the physical 
identification of solutions easier as we shall see in the 


next chapter. 
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The ultimate aim of the transformations is to 
describe the motion by a one-dimensional wave equation which 
can be solved for several velocity-depth structures. We 
introduce a new variable ZT that will measure depth in terms 


of density changes with depth: 





Zz 
rte f p(z) a. (2.14) 
: p 
O oO 


Si p@)) is a constant value (Phinney, 1970). We can see 


that for constant density f and (z-z ) are identical. The 
fe) 


equation (2.13) becomes 





ya s 
eT ee Bo(w) 6(t-2_) (2.15) 
dz p() 
Po 
q = (zy Q can be interpreted as a vertical wavenumber that 


includes the change of density. 

The homogeneous part of the equation (2.15) is the 
one-dimensional wave equation. If the velocity-depth 
dependence is such that an analytic solution can be found, 
the solution of the equation with the source is easy to 
find. The homogeneous equation has two independent solutions 
and every linear combination of these is a solution too. 

The equation (2.15) has two approximate WKBJ 


solutions 
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C 
A 1 ea if qdZ 
Pa i265) “= e (2216) 


I 
“2 


a 


except in regions of high velocity gradient, where q(Z) 


varies rapidly, or near the turning point f, , where q(Z_) = 0 
P T T 


wr 


(Morse and Feshbach, 1953; p. 1094). The factor q *iin ‘tthe 
WKBJ solution is such that the energy flux P * Vv remains 
constant in the vertical direction. Mathematically, the 
validity of the WKBJ solutions is explained in great detail 
an Budden (1961; W133)! 

The source lies in a region where parameters of 
the medium vary slowly and the two independent WKBJ solutions 
P, (5) and Po (5) represent waves travelling in positive and 
negative directions, respectively. We expect that the signal 
propagates along the same path from the source to the receiver 
as if it propagated from the receiver to the source. Hence 


the particular solution of (2.15) must be symmetrical with 


regnect tort and Ce ; Such a solution is 
BGreg ee Cr (2.17) 


where G is the lesser of (S55 .) and Gy is the greater of 
(o,0.)- This is a standard method used in electromagnetic 


field theory to find Green's functions (Jackson, 1965; p. 84). 
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The constant factor D is implied by (2.15): 





a7P. & 2. 
pe Ges oe Sia 5) p 3 6%) 


where € > 0 is an arbitrary real number. From the limiting 
case where € + 0 we obtain 
p 
O 
— FP (Ww 
5 () 


Ps 


to ~ ————— (tea) 
many Coates Ce) 

where W(P, (Zs Po (E)) is the Wronskian. In homogeneous 

medium, , and Ps are the exact solutions of the wave 

equation and the result (2.17) together with (2.18) is 

equivalent to solutions obtained by other methods (Ewing 

étwatasd19FAyvep lL ibd: «Chapnam and Phinney, (1972). Richards 

1970) has shown that in slightly homogeneous medium the 

Wronskian can be obtained from the approximate solutions valid 

in the source region. Thus the particular solution may be 


written as 


“w “a 


p “a “A 
oF , Wy SP, (GY oP, (c,) (2.19) 


oO 


P(E. ae oe oS 
2 P. 


The general problem of a source in inhomogeneous medium is 


a complicated one. Effects of both the point source and the 
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inhomogeneous medium exist together and complicate the 
solution. Ward (1973) used the principle of seismic recipro- 
city to solve a problem reciprocal to a problem of a point 
source near or within a linear transition layer. In the 
reciprocal problem the effects of a point source and the 
medium are separated and the solution is relatively easy 

to find. In the problems discussed here the source is 

always in the nearly homogeneous part of the model and we 


shall use (2.19). 


2.4 The response integral 


The source solution (2.19) represents waves 
travelling away from the source in both positive and negative 
¢ directions. We are going to study the waves reflected by 
an anomalous region below the source. Therefore only the 


downward travelling solution will be considered: 


pees Sak ag Mpaea etch (2.20) 
s 2 0 s uf C 2 ia ‘ 


Ss 





Its amplitude at the source is 
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While propagating downwards this wave changes by interaction 
with the medium. If reflected by the anomalous region, it 
propagates upwards towards the receiver oe? if transmitted, 

it continues to propagate downwards. Therefore, there are 

two waves above the transition, the incident and the reflected 
waves, while there is only the transmitted wave below the 
transition. The generalized reflection coefficient RCO 525, 5K) 
is a result of these boundary conditions and thus it describes 
the response of the medium. Hence the transformed pressure 


response is 


tJ >> 





; p ay. 
B(dytjoay =r5 steSiB qu) j PLCS) PACOs) RCCL.0 4K) Jo (kr)K dk 
12] 


(Jaen) 


phe result “of (2.22) is standard in seismology and has been 
used by many authors e.g. Ewing et al. (1957, ps 94) “and 
Chapman and Phinney (1972). The reflection coefficient 
R(G5o.>K) will be determined in the next chapter from the 
exact solutions of (2.15) for Epstein velocLty structures. 
Then the response integral (2.22) can be evaluated asympto- 
tically or numerically. Both methods are equally important 
in the theoretical investigation of spectral amplitudes. 

The first, while of limited range of validity, indicates 


the most efficient way of performing the numerical integra- 
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tion in the complex Kk plane. It also determines the correct 
asymptotic behaviour of the numerical results. The numerical 
method, while valid in the entire frequency range needed for 
computation of synthetic seismograms, demands sophisticated 
computers and is still rather costly. Richards (1973) 
Suggested a middle way between the two. His approach involves 
the numerical integration of (2.22) but the solutions used 
are the WKBJ solutions rather than exact wave solutions. Lt 
contains the higher order terms of the second order saddle 
point method but ignores the higher order terms of the wave 
solution. For the problems we shall consider the wave 
solutions are simple analytic functions so we need not resort 


to this approximation. 
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CHAPTER 3 SOLUTION OF THE WAVE EQUATION BY EPSTEIN'S METHOD 


In the previous chapter the equations of motion 
were established and reduced to the one-dimensional wave 
equation. We are going to investigate the wave propagation 
in media with continuous velocity transitions. The non- 
geometrical effects at caustics, critical points and shadow 
can be studied for the Epstein velotity atructures (Fit. 
3.1). We demonstrate here the method of circuit relations 
first used by Epstein (1930). His approach was discussed 
later in books by Budden (1961, p. 369), Brekhovskikh (1960, 
p- 172) and its application to seismic waves is due to 


Phinney (1970). 


3.1 Statement of the problem 


We will assume that the velocity-depth structure 


is one of those drawn in Fig. 3.1. The transition zone is 
centered around [ = 0 where the veloeity is Mee. Its thickness 
is measured by factor o so that § = - - The transition zone 


grades smoothly into homogeneous half-spaces above and below 
the anomalous region and velocity there approaches constant 


values v, and v 


l >? respectively. Thus we expect that in these 


regions the downwards and upwards travelling waves can be 
distinguished. This is an important property of our medium 


as it will help us identify the right solutions as waves 
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incident, reflected and transmitted through the anomalous 
region. 

With the new variable for depth 3 = ~ the equation 
to solve is 


“w~ 


cars o°4' (4)? “0 (3.1) 


This equation has two approximate WKBJ solutions (2.16) 

which represent two waves travelling in the opposite direc- 
tions and are valid in the regions above and below the 
anomalous region. We assume the source placed above the 
transition i.e. no wave travelling in the negative e direc- 
tion exists far below the transition (Fig. a 4¥, The solution 


that satisfies this boundary condition is 


: 19 (Aaa -i0 (4a 
2 ES ee for 0 >>4 


P(w,k,4) 
aye 


~ io daa 
P{W 5K 4) av oS T é d 
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ee 


and every exact solution of (3.1) must asymptotically approach 
to (3.2). Coefficients R and T are the reflection and trans- 


mission coefficients. It is therefore essential to find the 
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Figure 3.2 The WKBJ solutions in the medium with 
velocity transition. 
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relationship between the field in the region of incidence 
and the field in the region of transmission. Mathematically 
this is equivalent to determining the correct analytic 
continuation of these asymptotic solutions outside the 
range of their validity. 

Epstein's method shows how the solutions above 
and below the transition can be matched and the eennect 


reflection and transmission coefficients found. 


ee Compatibility of the wave and the hypergeometric equations 


Epstein (1930) showed that the hypergeometric 


equation 


2 
ey ee esa) cy 22. = ake <0 (3.3) 
2 d€ 
d€ , 
is equivalent to the wave equation 
2 
oa +ny=0 (3.4) 
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nh te (e-1)* f (c-1)*-(a-p)? *. (at+b-c-1) (a+b-c+1) ae 
4 (ed +1) 4 (e441) 
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Pe ee YE eer a ES ES (3.75 
. oi 6 (e441) 3 (e441) * 


The hypergeometric equation has three singularities, €.'=-0, 


1 
S5 mee 63 =o . \hhe\ohanva ‘ot variables (3.5) maps the 
inside of the unit circle in & plane on the left-hand half 
plane of 4 - The outside of the unit circle in & plane maps 
into the right hand half-plane aN Eo ssa) Bs Ena Oe 


The solutions ®, CL O22 354,556) are alf® in terms 


of hypergeometric series of type 


(a)_ (b) 
LS Se ee ae Sn = okaroy 
Pea bees e) = va cr n! a Ca). ad T(a) We BAT 


and are all defined in Appendix A.l. They are valid only 
in those regions of & plane where the series are convergent. 
If we map these into 4 plane we have the following solutions 
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For A] = 0 


¥3 (4) = Phe srnsh Fieve") e 2(c-1)2 (rety® 


¥, 9) : F(c-b,c-a;c-a-btl;lt+e?) e 2(c-lg (1+e4) Cab +4 
(33) 

For 3 Ee 

¥5 (4) . BOecet heise.) eo 2(a-b)4 (1te 4) 4 


¥6(g) = Fla,a-ctlsa-b+1;-e°4) eo #(- PY (14.74) 4 


Outside the region of convergence of the series, each of the 
solutions can be analytically continued into the rest of the 
g plane. The solution is then expressed as a linear combina- 
tion of the two series solution convergent in that region. 

The analytic continuation is possible because an integral 

solution of (3.3) exists that is valid throughout the € plane 
with exception of the singularities, ae Eo and 6. (Appendix 


A.1). The relationship between the solution in different 


regions of the convergence is given by their analytic continua- 


tion can be written in a simple matrix form: 


ler ar i # j (3.10) 
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12 is a 2 X 1 vector of two independent solutions in the 
region of convergence around oa: 144 Le a 2 X<2 matrix of 
analytic continuation from region by to region by (Appendix 
Aw.) 

We are looking for solutions of equation (3.4). 
These are connected to the solutions of (3.3) Dy (3.560 5 
Therefore the same matrix relationship as (3.10) is valid 
for solutions res 1 is a vector of two independent 
solutions in regions of convergence Sh aay Bet (41 a) 
32 = 0, Ba = +o). By inspection of the asymptotic solutions 


(3.2) we conclude we need the following relationship: 


% 74 PE bart 
& 2 (3.4151) 
x6 Sede etek *9 
The matrix is the continuation matrix 314 given by (A.1.11) 
Tr (1l-c)T(1+b-a) ['(c-1)0T (1+b-a) 
lr (l-a)T(1+b-c) r'(c-a)I(b) 
314 = es tt Ge 
['(1-c)T (1+a-b) l'(c-1)T (1+a-b) 
['(1-b) [(1+a-c) T(c-b)T (a) 
If the solution pha for 4 > 0 is given by ar or Pe and for 
3 < 0 by the linear combination of a and Py prescribed by 


the matrix of analytical continuation 314 we can be sure we 


Y and 


have chosen the right solution. Identification of Ya D 
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¥ 5 or ¥¢ as incident, reflected and transmitted waves depends 
on the choice of Im(c-1). If Im(c-1) > 0, the leading term 
in ¥, (3.9) gives the wave travelling in the positive Z 


direction and in our problem ¥, is an incident wave. ¥ 


iz 


represents the wave travelling in the negative 4 direction 


2 


i.e. the reflected wave. If Im(e-1) < 0 the directions of 

propagation of Ya and Fo would be interchanged and ¥, would 

represent incident wave while ay would be the reflected one. 

The sign of Im(a-b) which determines the wave propagation 

for ‘ > 0 has to be taken in agreement with the sign of 

Im(c-1). In case of constant velocity throughout the 

medium the wave travels without change. Therefore if Im(c-1) 

woeett 16 ¥. that represents the transmitted wave. ¥6 

represents the transmitted wave if Im(c-1) < 0 (Fig. 3.4a,b). 
The matrix of analytic continuation represents the 


law of reflection when the incident wave propagates from -~, 


The inverse matrix 


represents the law of reflection 


Rites 


Should the incident wave propagate from +“ (Fig / 3.:Acia). 
There is, however, a better way or more common way to write 
the reflection law once we choose the sign of IM(c-1). The 
equation (3.11) represents two physical problems that never 
arise at the same time . This is due to ambiguity in sign 
of Im(c-1). Once that is chosen, one of the equations in 


(3.11) has no physical meaning. The same is true should the 
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Figure 3.4 Analytic continuation interpreted physically as 
reflection-transmission law. 
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source be positioned at +”, 


For Im(c-l1) > 0 and Im(a-b) > 0 we can write 





eT Si4 a a) 
a5 a tse ia) 
if S54 :,,/ t. 
n 
52 1 
where § = — S = —- —— 
LTR a 12 es 
‘ detese Tea Ae 
21 V4 22 ro 
PUR ek? 
sigh’ edit mee 


The vector on the L.H.S. represents waves travelling towards 

the transition i.e. incident waves from above and below. The 
vector on the R.H.S. represents the waves travelling away 

from the transition after reflection or transmission (FPigx 3u5)% 
The diagonal terms of the matrix § are the reflection coeffi- 
cients, the off-diagonal terns are the transmission coefficients. 
We will show that suitable normalization of the vectors with 
respect to constant energy flux in piste Terie ays will yield 
identical transmission coefficients in both directions. Then 
the matrix S is symmetrical. We have found a new formulation 

of our problem which is analogous to that used commonly in 


seismology (see for example Chapman, 1971, unpublished lectures 









7 
a . 
‘7+ is bsnoliteog sd sa21u08 


me 


r « 
, A » «& s 
—-  « é — eS Tee 
7 ++ rf. 4- 
af 
Ps A ' ' 
é a ee 
: Be x 
r _ - = 4 
> &« . a re” 
+5 A i 
z 
7 
. 
dl = 
i 


4 @ingas sy ,2@.8.,d sf3 go toJ9sv & 
: : ; e. 
av mox? 3 raw Jneblant .s.!l notiteaat?s & 
oie 
| 
anilisvaz3 svysw oct etasesiqsz .2,838 ef3 ao aa 
a 


-— = 


eCC.& .gltt) aol tmene12 10 notisoslie:t redjis aolstenax? sd2 1 






$09 molijioesliery edi ste & zizijam $43 Yo ears? Kanoquad 9 





ztieo. Holeelauasyd of3 sis awse3 Lemogatb-21d oz 4 et 





35 


MEDIUM 1 





MEDIUM 2 


Figure 3.5 Reflection laws for incidence from above and 
below the interface. 
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for ‘Ehyetes .521).<cLt emus tcbeotséd -with care, however, as all 
the reflection and transmission coefficients are valid in 


far regions. The field 


ES ['(e-1)P(1-a) I (1+b-c) 

Sigh Whi bs Rhhee rca ET Coes mae) sas 
is valid throughout the region above the transition j1.¢. for 
4 < 0 without restriction. The field 

[(l-a) « ['(1+b-c) 


Oye tr Gite-ayraac) '5 3? hee. 


is also valid for 6 > 0 with no constraints. The reflection 


coefficient 


R( ) = '(c-1)T (1-a)T(1+b-c) : *2 Gr? es 
peter oi (e-ayithiTtize) 7G.) ; 


and the transmission coefficient 


P(l-a)T(itb-c) . *5 Fe) 


E Ba Ses PCiabda) Mie) ¥G,) Sete 


however, are valid only in regions where we can be sure that 
Pie ¥, and Pe represent waves incident, reflected and trans- 
mitted, respectively. And this is only in the region where 


medium is nearly homogeneous. Throughout this work we have 


found that the reflection coefficient and transmission 
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coefficient are accurate enough for frequency range and 
thicknesses considered. Comparison of the leading terms 
in the solution with the next terms gives the estimate 
of the error (Appendix B). 

To apply these results to equation (3.1) we must 
assign suitable physical meaning to constants a5. 05 (Ss 
Equations (3.1) and (3.4) are equivalent if Oq = h. Solutions 
(3.14) and (3.2) of the equations (344), and.(€3 <1), «must .be 


equivalent. By inspection of (3.2) and (3.14) we get 
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v +5 fe) 2 5 
a-b = 2ioq, = iS Bo foo o es: 6.) | = 240 — C - K*) 
2 fe) 2 : 2 
2 28 eee 
2 2 
3 ck 
where § is the effective thickness 
s- i+ 20-4. 26 
1 4 
From ong? = h*(0) we obtain 
+ ae Dy ed 2 Z *4 
atb-c = + 2y = + (1 + 160 (a5 - 6(q5 + q5))) (3,.18a) 


The sign is arbitrary and each gives the identical reflection 


and transmission coefficients (Appendix B). The analytic 
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expression for vertical wavenumber for the Epstein models 


sk 


2 
24 AP) 
2 2 ae 2 
q (a) = q, - (q)-q5) % 4d). - rem rds as Me 
d 4 e Shce = : FER 
(3.19) 


and its form is shown in Fig. 3.1 for various values of Vie 


V5 and we ; The law of reflection with normalization of 


solutions to preserve constant energy flux in yi direction is 


a= 


a5 ? " 2 
me 5 12 ay 
qd qy 
* = (3.20) 
“6 S S Me 
& “ps 22 “ay 
qo q9 


where matrix S is symmetrical with elements 
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If we use the asymptotic expressions for "469g? and ¥5 4) 


in (3.16) we obtain the following reflection coefficient: 


a Eee , oe Ie 


where Gy ei 


pire oe ae cy Lo 


) and 


The term in the exponent gives the change of the phase of the 


wave on its path from the source to the receiver. Using 


asymptotic expressions for ¥ 


1%»? and ¥5 (44) in (3.17) we 


get the transmission coefficient 


where 





-i10(q @ = q. 9 ) 
Wee Eee aa (3:23) 
Ll. 
(fe a & 
eee: RD 
Ve 


The reflection coefficient (3.22) will be used 


for evaluation of the pressure response (2.22) for the 


monotonic velocity transition (Chapter 4) and for the 


velocity reversal (Chapter 5). 
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CHAPTER 4 CAUSTICS AND HEAD WAVES DUE TO MONOTONIC 
VELOCITY TRANSITIONS 


4.1 The reflection coefficient 


For the monotonic velocity transition in Figures 1.2, 


Pasoand, 4, 1 a = 5(qy + a5) and the vertical wavenumber is 
given by 
2 oad ple 2s yt eee: ca 
q (cy = a(q, + q5) (qd) q,) tanh IG (4.1) 


The reflection coefficient (3.22) becomes 


3 2 : ‘ 
1,74, M(1+2iogq,) r (1-io(q,+q,)) a 8th 544, oy) 


Beer) ae, Fcdtioa,y 





r(1+i0(4q,-4,)) 


(qc. 


The reflection coefficient (4.2) was studied by Epstein (1930), 
Rawer (1939) and Phinney (1970), although they only studied 
its modulus and not the full wave response. The term 

P(1+2iog, )/P(1-210q,) effects only the phase change that 


occurs during the reflection. When attenuation is not 


considered gate) is real. For angles of emergence less than 
Vv 
the critical angle a = sin - 4 » the vertical wavenumber 
2 


below the transition is real and partial reflection occurs. 


In this case 
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1.6 km transition zone at 20 km 
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Figure 4.1 The partial reflections and head waves near 


the critical point for monotonic velocity 
increase. 
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sinh (wotqs-4q.)) 


Fdat gasa Gro (a, #45)) 


(4.3) 
At sub-critical angles the partial reflection depends on the 
erfective thickness S$ (3.18). The reflection and transmission 
through the medium depend on the absolute thickness of the 
layer 0 as well as the wavelength of the propagating wave and 
the effective thickness S includes both Slrects, “LF 3° Le 
small the reflection is nearly as strong as the sub-critical 
reflection for the sharp velocity increase from V4 to Vo: 

Then the reflection coefficient is the well-known Fresnel 
formula 

1g88 9 


[R| = 
44745 


(4.4) 


In this case the reflection coefficient is independent of 
frequency and the Cagniard-de Hoop method may be used for 
evaluation of the response integral (Cagniard, 1962; de Hoop, 
1960). 

If the effective thickness is large the nature of 
the partial reflection is changed. The partial reflection, 
which is quite strong at near-vertical angles when the transi- 
tion is thin becomes very weak for thick transitions with 
triplications (Figure 1.2). Most of the energy is transmitted 


through the transition and only in the neighbourhood of the 
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cusp some energy returns. The decay of the modulus of the 
reflection coefficient into the region before the critical 
point is given by 


le) v2 5 
2 

-"S5e— cet ~ ain O29 
a 

IR| =e (45> 


NSO bh 


For angles of emergence from the source exceeding the critical 
angle, |R| = 1 and total reflection occurs whether or not the 
transition zone is thin or thick. Thus the reflection coeffi- 
cient (4.2) gives the phase change along the wave path. 

The study of the analytic behaviour of the reflection 
coefficient is necessary for the evaluation of the response 
entegralwo2.22). Only perfect understanding of the reflection 
coefficient in the complex plane of the horizontal wavenumber 
K enables us to perform the evaluation of the response integral 
efficiently. 

The reflection coefficient (4.2) has eight branch 
points #K, = + , ke, = +2, 4x ee ond te mee BO 

2 s r 
due to zeros of the vertical wavenumbers dy* do» q, and qd, 
(318) and (3.22) 2 respectively. The complex k plane has, 
therefore, sixteen Riemann surfaces separated by branch cuts 
but only four are important in our discussion (Figures 4.2 


CV Oe oe » ee The twelve remaining sheets need never be considered. 


In our problem the best choice of branch cuts are lines Liars ; 
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A a ; Where p «= 1. 2, a or r. The 


branch cuts from Kie KS and ee to -K 


—“K , he respectively, 
are very close together and when speaking about the top sheet 
of qd, we implicitly assume also top sheet of a. and qb The 
bottom sheet of qi is also assumed to be a bottom sheet of qd. 
and : Oe The notation (++) sheet has been used by several 
authors including Phinney (1961) and Gilbert and Laster (1962). 
Their definition is slightly different from ours which suits 
our problem best (Appendix C) . The (++) indicates the signs 
of the imaginary parts of qo and qi> respectively. The (++) 
sheet is sometimes referred to as the "top" or "upper" Riemann 
sheet in this thesis. 

The reflection coefficient has poles at the poles 
of the gamma functions in the numerator of (he). - The 
contribution from poles of reflection coefficients often has 
important physical interpretation. For example: Lamb (1904) 
interpreted the residue at the pole of the reflection coeffi- 
cient in a homogeneous half-space as the Rayleigh wave. Pekeris 
(1948) interpreted normal modes, i.e. waves due to interference 
within layers of a layered medium, as contributions from poles 
of the reflection coefficients. Phinney (1961), Gilbert and 
Laster (1962) and Chapman (1972) have shown how poles on 
other Riemann sheets may cause separate arrivals or modify 


geometrical ray arrivals from the top sheet. We must therefore 
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investigate all the poles of the reflection coefficient. 
There are three characteristic or secular equations which 


determine the position of poles in our problem: 


2ioq, = -n RN toes Rte Se oe (4.6a) 
io(q,+q,) =n ft Letsas «ke (4.6b) 
q,+4, = 0 (4.6c) 


The first secular equation gives poles on the top 
Riemann sheet. They lie on the real axis "beyond" the branch 


Bownce Ik, (Pig. 42a) 


a 
a3 Ni + (> + 5 7) (HEF) 
vi 40 Q5 


The signals arising from their contribution decay exponentially 
in the vertical direction and propagate with phase velocity 


smaller than v.: 


1 
2 

2 p, —4s 
. ee Cee eS 
Se rorgh itite—5 Pm: 2) 
mad Py 


We are mainly interested in waves propagating with phase 
velocities between v4 and V5 i.e. arriving earlier than signals 


arising from the pole's contribution. 
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The secular equation (4.6b) gives sets of poles 
either on (-+) or on (--) sheet. We will treat here only 
the case of constant density as the algebra is simpler and 
overall conclusions are analogous. The poles on (-+), (+-) 
and (--) sheets must be investigated as they sometimes 
contribute signals called "leaking" modes. In contrast to 
signals arising from poles Ka? which decay with depth, signals 
arising from poles on other sheets grow exponentially with 
depth. The energy "leaks" into the lower medium. The 
"leaking" modes arise from poles that are within certain 
neighbourhood of a branch point and, therefore, influence 
the reflection coefficient on the top Riemann sheet (Chapman, 


1972). In our problem, the poles 


wi" te "e a" B- ae 2 *3 
oa L. 4 a See . SE fen 
Te Se {4 ( ig as mY + te Gs dant (4.8) 
v4 V5 fe} Vi V5 


would have to approach the branch points +K 4 to cause such 
effect. This is very unlikely as the poles lie always on 
the real axis "beyond" the branch points + Ky (Fig. AsBoe 


At low frequencies all poles lie on the (--) Riemann sheet. 
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As frequency increases the poles move towards the branch 
points and a cutoff frequency exists for each order n when 


+ KO appear at the branch points +kK, and pass from (--) on 


2 i 


(-+) sheet: 


gee | ui 
Bg ag 9) 
sO Te 
For higher frequencies the poles t5KO move away from the 


branch points on the (-+) sheet (Fig. 4.3a), Thus the two 
sets of poles tok will not become important when evaluating 
the response integral (2.22). 

The roots of equation (4.6c) give poles analogous 
to those in two fluid half-spaces, Gilbert (1964). They lie 
on the real or imaginary axis of (-+) and (+-) sheets: 


22 
04%o 








(4.9) 


Depending on the parameters Py» Po and v the poles can 


1? Vo» 
contribute to signals propagating along the thin transition. 
We are interested in the propagation at the receiver point 


well above the transition and the contribution of these pulses 


is negligable there. 
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We can conclude the investigation of poles by 
saying that for the monotonic velocity transition, the 
signals arising from poles" contribution arrive at later 


times than those we are going to study. 


4.2 The WKBJ reflection coefficient 


In Chapter 3 the reflection coefficient for Epstein 
models was established by solving the wave equation (3.1) 
exactly. The approximate WKBJ reflection coefficient can be 
obtained from the WKBJ solutions (3.2) and it is given by the 


well known "phase-integral" formula 


R(T +o, 5K) ah oe = (4.10) 


(Budden, 1961; p. 329). The integral in the exponent may 
be regarded as a line integral in the complex Zt plane around 


the branch cut given by Im q(t) = 0 (Fig. 4.4). is the 


or 
turning point of the ray and the formula (4.10) can be deduced 
by solving the wave equation in its neighbourhood. The WKBJ 


solutions break down in this region and must be replaced by 


solutions valid there. For linear behaviour of of 
) C4, 11) 


the wave propagation in the vicinity of the turning point is 
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described by Stokes equation 


7 J i 
ot - xP = 0 (4.12) 
dX 
243 
where X = fa (tT - Cr) r 


Solutions to this equation are any independent combination 
271i _2Ti 
of Airy functions Ai(x), Bi(x), Ai(xe 3 ) or Ai(xe 3 ri (See 








Appendix A.2). The WKBJ reflection coefficient (4.10) is 


implied by the formula (A.2.2) if we assume |x | large. Then 


277i Prue! 








Ai(xe ° ) and Ai(ye ? ) represent waves travelling downwards 
and upwards in regions above the turning point. In the region 
below the turning point Ai(x) describes the evanescent wave. 

The WKBJ solutions for real and imaginary wave- 
number q represent either oscillatory or exponential solutions. 
They remain valid for complex wavenumbers (Budden, 1961; Pr 435736 
This will occur in the response integral (2.22) when the 
stationary phase point corresponding to the ray is off the 


real axis in the complex kK plane. The vertical wavenumber 


2 
q = Ay - “2% cannot have zero for complex K unless the 
v 
velocity v = v(t) is analytically continued to the complex 


GC plane. Then a complex depth On exists so that q(on) sate Sa 
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lines of the Airy functions, (see Appendix A.2), are moved 


in the complex ¢f plane. If we define 


103) 


Qa = q dz (4.43) 


the Stokes lines are defined by Re a = O and the anti-Stokes 
lines by Im a= 0. For complex wavenumber kK they start at 
complex depth Cp q(.) =i tt In @.> -0- then Tm or m1) 

as follows from the Cauchy-Riemann conditions (figs 4a b)x 
Both the downward and upward travelling solutions are 
dominant on the real ¢ axis (Fig. 4.6b), i.e. exponentially 
growing in their direction of propagation. If Im « < 0 then 
Im Cr > 0 and both solutions are now subdominant on the real 
G¢ axis, and, therefore, they both decay in their direction 

of propagation (Fig. 4.6a). This case also corresponds to 
the physical problem of absorption when the wave is partly 
attenuated in the medium. Ewing et al. (1957, p. 272) showed 
that attenuation can be theoretically expressed using complex 
wave velocities. Thus for real wavenumber kK the Fig. 4.56 
applies to this case. If the wavenumber k is complex, the 
presence of attenuation reduces Im Cr . (Hi flertdively, it 
decreases the exponential growth of WKBJ solutions if Im kK > 0 
and increases their exponential decay if Im K < 0 (Heading, 


1962: pp, 75). 
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The complex f plane for K complex 
(a) Kem im O,. Ch) Intth Sok. Iel he 
Stokes (S) and anti-Stokes (A) lines 
are also indicated. 
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{a) 
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UPGOING = SUBDOMINANT 


(b) 


DOW NGOING = DOMINANT 


UPGOING=SUBDOMINANT 


Figure 4.6 The Stokes diagrams describing the 
behaviour of the WKBJ solutions. (a) 
foraiasts. 2° 0. ib} yeger In &y < OG. . When 
inside the circle the heavy line indicates 
the subdominant solution and when outside 
the dominant. 
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The approximate reflection coefficient for the 


monotonic velocity increase can be derived from the phase- 


integral result (4.10). The vertical wavenumber 
+5 
2 2 2 2 G 
= 1 Sa os a 
q(T) (aCq, + q,) a(q) q5) tanh ao? (4.14) 
has branch points at q = 0 and q = © 
2 , 
sad! Sis 1 
= t A P ale 
4s 2o {arctanh ( 7 >) + i(2j + 1) a 
Scieatay: 
tC. = i+ 20(4 +1) 5 j = 0,41,+2 
co? Fj J 2 ’ I s—-+s—-“s eee 


ee bs 


respectively (Fig. 4.7, Budden, 1961; p. 448). 

From analogy with the linear behaviour of ante) 
we expect that one of the branch points aha? = 0 above the 
real Cf axis must be considered for evaluation of the phase 
integral. We perform the integration along the branch cut 
Im q = O starting at the source Ce? where q(t.) cde Ui + ES 
and ending at the receiver ce on the other Riemann sheet where 
a(S) me: te e.. - This is represented by curve Ci in the 
complex € plane or complex q plane (Fig. 4.7). In the 


complex q plane the curve Cy can be moved to our convenience, 


provided we take into account all singularities enclosed or 
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COMPLEX ¢ 
PLANE 


COMPLEX aq 
Req PLANE 


The complex ft and q planes showing the 
contours of integration used in the 
evaluation of the WKBJ reflection 
coefficient for Epstein profile (Budden, 
1961; p. 448 and 449). 


ee 





> kaJIMOD 
IAAIS 


op AZISMOD ’ 
(IL 
IVA S49 


60 


crossed by the deformation of the contour. The best integra- 
! 


tion path in the q plane is C, because the integration is 





1 
then very simple. For Cs = Cy we get 
q,7€ 
: 2 2 
Ly q 
q dt = 20 ( 5 5 + 5} x) dq + 
es 58 roks 
Cy €-q, 
49 
+ 220i Res ¢ =A ) 
2 q=-4, 


and finally 


(4.16) 


This result obtained from the phase integral is identical to 
the approximation of the exact reflection coefficient (4.2) 
for large effective thickness. Then we may use the Stirling 


approximation for gamma functions (C.11) and 


Piliectie se) Pee te easae (4.17) 
where O(K) = ; q dc for Ce = ce . 
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The phase $(k) is real for angles of emergence from the source 
greater than critical. The phase change > is in agreement 
with the phase change due to total reflection.© The* two 
logarithmic terms can be interpreted as phase change due to 
the transition layer and 2415, is the usual phase change along 
the wave path. 

The fact that the WKBJ reflection coefficient agrees 
with the high frequency approximation of the exact reflection 
coefficient justifies the choice of the branch point a for 


evaluation of the phase integral. The branch point of for 


1? 
example, would give a wrong constant of integration due to 
encircling both poles q = qd5 and q = “45 (see curve C, an. 
and q planes in Fig. 4.7). The singularities must be well 
separated so that they do not influence the evaluation of the 
integral. This requirement, however, is implicitly satisfied 
as the WKBJ solutions are valid in slowly varying media i.e. 
in our case broad transitions. More detailed description of 


phase-integral method may be found in Budden (1961, ps 437) 


and Heading (1962). 


4.3 Asymptotic behaviour of the response intepral 
rn ee VE response integral 


For the monotonic velocity transition the response 
integral (2.22) is 


+00 
; re) 
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P(w,r,t_) i c P .(w) R(T 505K) Boo? (e¥) 1, dk (4.18) 
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where R(o55.>K) is given by (4.2) (see Appencix C). The 
integral (4.18) can be evaluated approximately at high or 
low frequencies. These approximations will show the best 
way of performing the numerical contour integration in the 
K plane in the entire frequency range. 

At high frequencies the WKBJ reflection coefficient 
from the previous section may be used. If we use the 
approximate formula for BS) (cr) at large arguments (C.4) the 


response integral becomes 








e +00 op or 
P (w) 0 -i- ikr-i qdt +i qdZ& hs 
no.r,c) = * a : e sae ae 
s+ = Ses 
t av2nr Ps Ss or Ts 
—©o 
T -” 
P (w) 0 -i- : bs 
~ _s 5, SEO ie 4 oikr + 10(K}) K- ae 
2V2Tr a q. 
= 00 
(4.19) 
where 
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Eas Siig 
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The integrand in (4.19) is oscillatory along the real kK axis 

and contribution from the integrand is significant only at 
points of stationary phase. The contour of integration along 
the real axis can be moved to pass through the stationary 

phase points (Fig. 4.8) and the integral evaluated approximately 
by the second order saddle point method (Morse and Feshbach, 


1953; p. 440). The phase in (4.19) is 


ny Gr 
f(K) = Kr - q dt | q dq = Kr + O(K) C4 S27) 
e or 
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Figure 4.8 Saddle points in the complex k plane for a_ broad 
velocity transition. 
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gives at least one saddle point on the real axis between K» 
and Ki: Cr is the depth of the turning point of the ray 
represented by the saddle point and it is a function of k. 
This, however, does not effect the result of the differentia- 
tion in (4.21) as q(,) = Up. silnecresult,in €4.22) can be 


compared to the classical result given by Bullen (1963, p. 111). 


We just have to transform his result to flat geometry to get 


horizontal range X(p) = r(K) of the ray that bottomed at 
depth Zp 
Zz Zz 
Ss r 
dz dz 
X(p) = p eres + p Soak 88 : (4524) 
en E-P es) te =p. 3 
map ay 
sin 6 ar? ? 
- it oe the ray parameter or horizontal slowness 
n= = is the slowness of the acoustic waves 


For the monotonic velocity transition where the triplication 
occurs (Fig. 1.2), three saddle points exist on the real 
K axis at ranges between O and 7 (Fig. 4.8). The saddle 


> W : 
point close to 7. represents the ray close to the critical, 








as okoq gatnzas od Yo doqeb o€9 BAD. 
fi 2 i bos taioq sibbsee ofa vd beineesiqges 


i? Yo sini $73 Jost39 tom 2aaob , tevaved 


> a = r 7” 4 a re a | 
: : ef 0 (25) ee (I¢.a3 

4 a - 

~ A esis si o3 E 
s i @ } a7r3 .03 ved 


, wole ! josiios xt9 aseJometag yer Sadi at 


Uiim 


eovew stieuecs af3 to eesnwote ad} et 






ramp — | 


soe. 






oO 





66 


the saddle point close to —— is the shallow ray and the 
saddle point in the middle _ the ray reflected by the 
high velocity gradient. At ranges 9 > O, only the near- 
critical ray exists i.e. the saddle point close to — : 
At ranges 0 < OG only the shallow ray with a saddle eee 
close to — exists. 

“a 

For narrow transition in Figure 4.1 the stationary 
phase condition gives two real saddle points at horizontal 
distances greater than the critical Xo (Fig. 4.9). They are 
the head wave and the reflected wave shown in Figure 4.1, The 
head wave is represented by a narrow saddle point that 
approaches the branch point — as the transition gets thinner. 
Concurrently, the steepest descent path through this saddle 
point approaches the branch cut integral which is interpreted 
as a head wave in case of the sharp velocity increase (Cagniard, 
1962). 

The stationary phase condition (4.22) is a trans- 
cendental equation and must be solved numerically. The phase 


at the saddle point kK* gives the arrival time of the ray 


C 


s o 
£(kK*) = q(K*)dZ +| a{e*) dt + Kee (ek) = wT (4.25) 
Sop So 

Again, T agrees with the classical result given by Bullen 


(1963, 0,132) 
















eda bee ver wolisdsa saa et ~ 02 seols gnkoqg olbbse 8a9 
I ~ 7 
oft yd bersetIte+: yar s#i of al bbim sad gl Jaieg stbts 


- 
-ysen 942 ylac a” < 0 sonnet 3A .doekberg vdisolsv Wyle 


— o3 seols Intoq slbb|es of% .&,.2 Seles Yar leolaths 

-* ork 
inicg sibbae «6 atiw yas wollede efa yiao ah > © segnet 2 
oe 0 


a 


-etetxe — 69° 98S 
- . 


rrecotiase of3 I.e siuatY at adlifenaxr? voxrvan yot ~~ 
_ 
[astnosiion Ja eintoq sibbsaa Isst ows asvlg notsibao0 oeadq 


os 
load 
¢< 
Q 


~(2.0 «tT 5* {eaitattip off aed? te3g9%R BSORST tb 


efiT .1.4 esuett at awode oVew beiseltay sh bas ovaw beod id 
( ; 


teadd tnloq sibbes worsen o ¥d batnaverqs? sk svaw Seer 


, b 
tesnnids si6g9 mwoltieandays 43 @s —_ jateg donaad sat asdosot 
c 


7 
albbae aids dguords daeq tnacesb Jasqes3a 8a .(isnerI94a0 
besterqisint ef dotdw Lleyvgeint jus sdoas1d oA3 sedosoxzgqs Ialog 


~~ 


»btelags)) sesstoot zitoolsv grade efi to 9889 at SVRW DEOR BG 
.” t cor 
(82 

ae 


-ens12 ® at (SS.4) aolstbno> seailg yraaoliase/sdT \ i? 












saedg sit «.yllaotismua beviow ot Save bas poliaups tngnebasi 


- - a lke yi a 
|. © gmx edd Bo eats Laydt2e ad3 covig #3 aekoq ofbbas, +4 
7 ee Sa ee ee ey 

| : ! - r. 7 ; en cake 





67 





Rex 


=e 





exe 





Ww 


— oe Rex 
Va \ ] 


ny 
7 


Figure 4.9 Saddle points in the complex kK plane for a 
narrow velocity transition. 





The second derivative at the saddle point is 
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dz (4.26) 


(4.27) 


(4.28) 


determines the direction of the saddle 


eee | 
* gives the amplitude of the ray. 


The forward 


branches of the triplication AB and CD are due to rays that 


travel further in horizontal direction if they penetrate 


deeper into the medium, 


The saddle points are positioned at 


Bd 
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T 
angle — with respect to the real kK axis since oa “COand For 
4 dk 


these rays the second order saddle point evaluation gives 


Pw) le) 4 =6iWT 
P(w,r,o_) =——-+-—- re (4.29) 


/ fe) 
evat g qd. {x (e) |< ig 





Thus there is no phase shift of the resulting response 

with respect to the source Es {t)} and the eur tess is 
frequency independent. The factor {x |= aS e represents the 
geometrical spreading of the ray tube along the ray path. 

The rays of the reversed branch of the triplication 
arrive at shorter ranges if they penetrate the medium deeper. 
The same is true for the reflected branch in Figure 4.1. 
Thus the saddle point is situated at angle - = towards the 


4 


: : dr 
real K axis since a > 0 . For these rays the pressure 


response is 


B(w,r,t) = Y icaaere re aes ial daa (he 50% 
q, fr («) | = 


and they are shifted by - 5 with respect to the source 
function P(t) - As frequency decreases the "geometrical" 


amplitudes (4.29) and (4.30) become poor approximations and 


(4.18) must be evaluated numerically. 
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The triplication has two cusps at B and C and the 


gurveie vs. 8 Ci.e, a vs. 0 curve) has two reversals at 
ot _ = 0 (Fig. 4.10a). The travel time curve for narrow 
K 


transition has a cusp at C which approaches the critical point 
of the sharp discontinuity as the transition becomes thinner. 


The *« ws. X curve in Figure 4.10b has»also a reversal at 


Ss 


Ac res - At all these points the second derivative of 
K 


the phase (4.27) is zero and the second order saddle point 
method cannot be used. The saddle points approach each other 


at ranges close to Oa» oO. or XG and separate evaluation of 


their contribution is not valid. However, the width of the 
anne dey he=. 35 
saddle point is proportional to |—= ie Vu) and, therefore, 
dk 
at high frequency the separate evaluation remains valid at 


distances closer to r (Kk) =r (r applies to either of On» 0, 
and X,)- The sign of the third derivative of the phase deter- 


mines at which endpoint the reversal occurs: 


<0 at the "near" endpoint 
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Figure 4.10 Horizontal wavenumber versus horizontal 
distance (a) for broad transition with 
triplication, (b) for narrow transition 
with head wave. The dots indicate the 
partial reflection. 
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When approaching the range ry from the illuminated region 
two very similar waves arrive. They differ only 
slightly in their path-lengths and amplitudes, they 
interfere and cannot be separated. This is caused by focusing 
of rays into a caustic due to the change in velocity gradient. 
Other examples of caustics exist in seismology which are 
formed by focusing of rays due to geometry as, for example, 
the PKP caustic in the Earth's core. At the caustic the 
saddle points coalesce, and at ranges out of the illuminated 
region they separate again and move to complex positions in 
the K plane. The steepest descent path passes now through 
one of these saddle points only (Figures 4.8 and 4.9) and 
its contribution gives partial reflection decaying rapidly 
away from the illuminated region. The stationary phase 
condition can still be solved for complex K and the arrival 
time determined. The geometrical arrival time is essential 
for the synthesis of pulses discussed in the next section. 
The response integral in the region of caustic can 

be evaluated by the third order saddle point method (Airy, 
1838). The first term in the asymptotic series was used in 
seismology by Jeffreys (1939): 

+00 


‘re ee 
PED aK = geet t )' gS £) y  aacyy) (4.33) 
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where 
= — e Le ee ae 
K | <a 
The definition of the Airy function Ai(x) is in Appendix A.2. 
The integral expression for Airy function (A.2.7) has two 
stationary points in region "before" the caustic (x < 0) 


which corresponds to the illuminated region: 


= tt fy y real 


Vy corresponds to the minimum of the phase i.e. to rays on 
branches AB or CD. Y» represents the maximum phase saddle 
point i.e. rays of the reversed branch BC. The two saddle 
points give rays that arrive to the receiver in the illuminated 
region at different times. No real stationary point exists 
in regions "beyond" the caustic (x > 0) where no geometrical 
rays arrive. On the caustic (x = 0) the two saddle points 
Ya and V5 coalesce at the origin. 

The phases at the saddle points before the caustic 
are 


£ (kc) ~ ly [3/2 


Fh 
F iia 
< 
loom 

Ve 
! 

Wit 


(4.34) 
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We have obtained an expression for absolute value of the 
argument of the Airy function in (4.33). yx = 0 corresponds 
to one of the geometrical endpoints of the triplication. 
Therefore if we know the phases of the two rays arriving at 
distance r(K*) in the illuminated region, we can calculate 
the integral at distance r' = r(k) “ r(k) me THUGS iw ohlie 
range r' is symmetrical to r(k*) with respect to r(k) . We 
denote T, and T, the geometrical arrival times at r(KQ) and 


C B 


r(ky), respectively. Ti> T, and T, are arrival times 
corresponding to branches AB, BC and CD at range r(x*), 


respectively. The asymptotic expression for the pressure 


nesponse near the caustic C.is 





2 Cr ee pra s Cc 4 MMe 
Memnxiee Figen Divig: 71) and ¢ 
4 2/3 
: Ai(tly wT, if T,)] ) (4.35) 


T, and T are arrival times at Rony where ry is the greater 





Of oa; ro + (r, - r)). Near the caustic B the response is 
T 
p 1 8 2 -1/3 iwT, - i- 
T. % d B 4 
B(w sr, 5.) msigte) go 365) Aue dé berpbheed of 
Ss qd. dk 


* Ai(t [+ w(T, - T,)] ) (4.36) 
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T, and T, are arrival times at me“ where re is the lesser 


gee(z. tT. + (r, - r)) . The - and + signs correspond to 


response before and beyond the caustic, respectively. ‘And 
2 
ee 


=| te 'given by C4s52) 4. 
dk 





The shift of the amplitude maximum from the geometri- 
cal endpoint r(kK) towards the illuminated region is shown 
in Figure 4.11. The peak of the constructive interference 
occurs before the caustic and causes the maximum amplitude. 
Figure 4,11 also shows the amplitude curves obtained from 
contour integration in the complex k plane. The agreement 
of the "exact" amplitude at w = 20 with the asymptotic 
amplitude (4.35) is very good. The exact amplitudes oscillate 
due to interference of the rays on branches BC and CD. Their 
geometrical amplitudes given by (4.30) and (4.29) are also 
indicated. The rays of the reversed branch BC are stronger 
and the amplitude curves oscillate around this value. Out 
from the illuminated region the amplitudes decrease rather 
quickly at high frequency. The contribution increases as 
frequency decreases and confirms the fact that partial 
reflection is enhanced at lower frequencies. 

Results (4.35) and (4.36) based on the Jeffreys' 
evaluation of the first term in the asymptotic series are 
correct when the two interfering rays are identical or, in 


other words, the two approaching saddle points give identical 
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contribution when considered separately. Which means that 
the K vs. X or O curves are Symmetrical around the reversals 
(Figure 4.10). This is the case of the broad transitions but 
as the transition becomes thinner the curve is distorted tit 
it coincides with curves for sharp discontinuity. We expect 
this to happen as it indicates the transition from the caustic 
behaviour to the critical point behaviour. The departure 
from the Airy phase behaviour can be compared in Figures 4.12 
and 4.13 which show the amplitude curves for two narrow 
transition thicknesses. For w = 110 the agreement is good 
for 1.6 km transition while at the same frequency for .2 km 
transition, the decay of the exact amplitude is slower. The 
partial reflection is stronger for thinner transition and 
results in enrichment of the amplitudes before the critical 
angle, In Figure 4.12 the high frequency behaviour given by 
the reflection coefficient (4.5) is indicated. It does not 
include the remaining terms in the integrand but its shape 
agrees well with the decay of amplitude at w = 110. 

The higher order corrections to the Jeffreys' result 
(4.32) can be obtained from the further terms in the asympto- 
tic series. These include differences between the interfering 
rays and were given by Chester et al. (1957) and Yanovskaya 


(1966). 
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The treatment of the caustic region described here 
is distinctly different from treatment of the shadow at the 
core of the Earth (Chapman and Phinney, 1972). The similarity 
of both problems suggested by the fact that two rays arrive 
simultaneously at the edge of the shadow, and at the far end 
of the triplication, is only superficial. In case of shadow 
the direct and reflected rays have different characteristics. 
The contribution from the direct ray remains finite as the 
shadow edge is approached while the reflected signal tends 
EO ¢ZCrn0. 

In spite of the limitation of the previous results 
to high frequencies they are very useful. They show that the 
easiest way to perform the numerical integration in the complex 
K plane is by integration along the steepest descent path. 
They also give the geometrical arrival time indispensible for 
computation of synthetic seismograms. The geometrical ampli- 
tudes are also used to form part of the spectrum at the high 
frequency end. This is discussed in more detail in the next 
section. 

At low frequencies a different approximation can 
be obtained. In section 4.1 we have shown that for thin 
transitions the reflection coefficient approaches the reflec- 
tion coefficient for a discontinuous velocity increase (4.4). 


This is independent of frequency and in such cases the 
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Cagniard-de Hoop method of evaluation of the pressure response 
in the time domain is possible (Cagniard, 1962; de Hoop, 1960). 
We shall investigate to what extent the method is applicable 
to our problem. 

The gamma functions in the reflection coefficient 


(4.2) can be expanded in power series according to 


1 a k 
eeteetie, FTAfnerx |x| < @ Lavoe 
l(1+x) k=0 k 
x is complex and a> Ley a, = (oe 2le } ay = -.655878 (see 


Abramowitz and Stegun (1965, p. 256)). This expansion gives 

a product of three power series in the numerator and deno- 
minator of the reflection coefficient. For |x| << 1 the 
series in the denominator can be expanded binomially. This 
restricts the validity of the following results to small 
thicknesses oO and low frequencies, i.e. to small values of 

the effective thickness S (3.18). In this case the reflection 


coefficient is 


q,24 Fig. (6 th) 
Bc ete a Tg gr oe 3 


(4.38) 
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We assume o << Poets, | and la, | << ph; lq, | Seibberethe (first 
term reduces to the reflection coefficient between two fluid 


half-spaces with velocities and densities v and v 


hnpeaaa a Pes 
For use of Cagniard-de Hoop method the results have to be 
reformulated in terms of Laplace and Laplace-Bessel transforms 


(2.9) and (2.10). We-.wére using the Fourier and the Hankel 


transformations as the inverse Fourier transform is easier 


to evaluate numerically. After the change of variables 
4 
1 ; 
-s = iw sp = ik -sn =- s(—> -p ) =f me q 
Vv Po 


the pressure response becomes 


b fas +i© 
Bete ie) ee sR ery f R(z_,2_38p) 0K (spr) & dp 
a aa T 2 Bier) re) ue 
-jo 
(4.39) 
where‘ for z = zg 
r Ss 
Map ariae p- ~“28T oC mw) 
Ee 2o a 2%) | aes 2 re) ‘Oils Pt 
R(z 5z_,8p) = ————— (1 - 4s“a (2a,-a,) —— * nu ina)e 
ol ae NyPotNoPy 7 wey f P4Po Leg 
=-2STl, (2E-27) 
2 
= RU(1 + £,(p)s7e Ee (4.40) 


The inversion into the time domain of the response integral 


(4.39) can be accomplished by method introduced by de Hoop (1960) 
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and used by Helmberger (1968). The paper by Helmberger is 
the best reference for details of the method used to find 
the results presented here. If we assume a Heaviside step 
function to be the source, i.e. P(t) = H(t), the time pressure 


response is 


Vv 
as Tor. 6.<, 64 = sine _i 
c v 
2 
42 
P(t,r,z_) = Pa(t,r,z_) + a (Pip (tsrs2_)) (4.41) 


ia for § > La 


P(t,r,z_) = Py fteots 2) + PR(ty»r,z,) + 
az 
+ 2 {Py y(t eotsz.) + Pig(ty>r.z.) 


(4.42) 


where t. is the lesser of (t=) and t, is the greater of 


Pee and R = ‘fm + 4|z Zz 125% 
Vv Bid 
a4 
E 
Dp R 
Pp(t,r,z) = 4 EN apes casein a 5) 
1 
‘ (t-T + 2pr) (tz) 2 (17-2 ) 
R 2 
ee 
t 
Py(tsr.z,) = 2 f psin 2p at (4.43) 
H r TT i, L oR 2 
ty (t-tT+2pr) “(t-T) (—5-T ) 
Vv 
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In, le 
where 8 = arctan (- ORES 
NyPo 
mer tees pee | (4.44) 


defines the so-called de Hoop contour in the p plane along 


which Imt = 0 (Figure 4.14). PR(t,r,z_) represents the wave 
reflected from sharp velocity increase arrivingjat t, = a 
1 


Pi (tsr,z.) is the head wave arising at critical incidence and 


propagating along the interface. Its arrival time is 








1 
r 1 bike 
cae v5 es he) 2|z-2. (4.45) 
sre ws 
The remaining terms in (4.40) and (4.41) 
a2 
Pi (t,r,z_) or (Pj p(tsr,z,)) for 6 < G. 
dt 
a? 
Pi (t,r,z_) = ne 1Po gt rer b se) + Rig (ty>r>z,) Lox 9 > or, 
(4.46) 
where 
p2 ‘ nas “Wis 
d 
Rip 2 = . 4(a{-2a,)o° - ; . j Re{ ae ay = 5 
= ee 2 
dai < (t-1t + 2Zpr) (t-1) “(1 -A5) 
R 
bas: 
Le 
Paes Te ee Po f nylnz ip cos 28 dt 
1h. x i oa. 0.2 . 5 1. Re 28 
H 2 
a 


(4.47) 
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are the corrections to the reflected wave and the head wave. 
They are due to finite thickness and continuity of the 
transition. In expressions (4.43) and (4.47) the parameter 
p i8 a function of Tt (Fig. 4.14), and all the integrals are 
just temporal convolutions that can be evaluated numerically. 
For computation of theoretical seismograms, this however is 
uneconomical, as the frequency range where the pulses are 
valid is very limited. The evaluation of the response in 
the frequency domain and the inverse Fourier transform are 
more practical and have wider ranges of applicability in 

our problem. 

The amplitude for a head wave propagating along a 
velocity discontinuity is inversely proportional to frequency. 
Thus even if the velocity increase is continuous we expect 
the amplitude behaviour of similar nature. For this purpose 
the amplitudes of the head wave were evaluated by contour 
integration and results are shown in Figures 4.15 and 4.16. 
We must distinguish two separate effects. One is the 
interference of the head wave with the reflected wave which 
changes the frequency dependence of the head wave within the 
interference zone (Cerveny, 1962). This éxphetad why the 
amplitude-frequency dependence approaches to that of “ as 
epicentral distance increases. The other effect is that of 


transition thickness. For thinner transitions the head wave 
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is more like the "real" head wave at the velocity discontinuity 
and, therefore, its amplitude becomes inversely proportional 
to frequency as the transition thickness decreases. At high 
frequencies the amplitude approaches the geometrical amplitude 
(4.29) which is independent of frequency. This is another 
phenomenon that does not exist at sharp discontinuities while 
the amplitude is inversely proportional to frequency in the 
entire frequency range. Thus, whenever the velocity change 
has finite dimensions, a finite frequency exists at which the 
amplitude behaviour changes continuously from = £01). ) this 
frequency depends on the dimensions of the transition and 

with decreasing thickness it moves to higher frequencies. 

This is also well demonstrated in Figures 4.15 and 4.16. 
Equivalent results were obtained by Lang and Shmoys (1968) 

who studied analytically high frequency asymptotic approxima- 
tion for the reflected field. Their approximation consists 

of two terms: the first is interpreted as a reflected wave 

for arbitrary transition thickness while the second represents 
a head wave when the wavelength is much greater than the 
transition thickness 0, and a reflected wave when the wave- 


length is much smaller than the transition thickness o . 
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4.4 Synthetic seismograms 


The results of the previous sections are not correct 
in application to wave propagation within the Earth if the 
transition is in great depth. Then the Earth's curvature 
cannot be neglected and the problem should be solved for 
spherical geometry. The earth-flattening transformation has 
been used by many authors (Muller, 1971; Hill, 1972; 
Helmberger, 1972) to obtain approximate solutions in flat 
geometry. The following transformation was used by Hill (1972) 


and applied to our problem: 


R 
z= Ro log rR 
r 


R 
v-(z) = cae v,(R) (4.48) 


R 
where (R,@,6) are the spherical coordinates and R. is some 
reference radius, in our case the radius of the receiver. 
v, (R) and v, (2) denote the velocities in spherical and flat 
geometry, respectively. Chapman (1973) showed that for a 
fluid the optimum transformation requires the following 


density transformation 


R 
p-(2) =e pb (R) (4.48a) 
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where p, (R) and OP, (2) are the densities in spherical and flat 
geometry, respectively. The earth-flattening transformation 
does not change the kinematic properties in the transformed 
model. The rays emerging at the same angle arrive at the 

same horizontal range at the same time. Thus the increased 
velocity in the flat model compensates for the longer ray 

path due to the depth transformation (Fig. 4.17). The dynamic 
properties, i.e. amplitudes do change and Muller (1971) showed 
that the ratios of the geometrical amplitude (4.29) to its 


equivalent in spherical geometry is 


sin 0 "6 
(A) 


oT > 
Hh 


10) 


Thus the amplitudes in the flat half-space are systematically . 
smaller than those in a sphere. 

Three models for monotonic velocity transition were 
chosen for calculation of synthetic seismograms. 


Model x of broad velocity transition in fig. Lvs 


ee 8.4 km/sec eee v , (6200) = §.4 km/s 
eee 10.4 km/sec Mr v (5810) = 9.74 km/s 
o = 35 
z = 180 km 
° 
z = z2 = 0 km R = 6200 km/s 
r s r 


p = constant 
(4.49) 


. sesh tants 


Slisiieve to notvalustas 169 nme 


h) ¢ (oose>, 


cuteay’ 


faolj0rom 103 elebom se tdT 









~~ 
_ 
1 e429 a78 (3).6 ban 
: 
“aqz7s edT -¥ievis; eqes? 






'GO7Iq 2Lismeactd eoff3 seoaedo jou 
oo3 38 gotgisae syat sdT 
an 7h 7 } a4 S2ns 7 fasjnzosizend s ‘ 
a, 
a7 ' . is _o 9nit? ty d3 ot vyitoo Ss 
‘ 
ba " a 
Ofjemioisgeets tiqeb silt of sub 









afs ob esbutilems .s.t , 8ottisgor 


(B82 TIOMOSS SHI Pa eolrazx sdZ 


= 4 
- 5 - °° 
7: 
a —_ 


sage @ mf saod? gaad3 x98 oa 






sosd to x lobow \ue 






/ 


Figure 4.17 


92 


Effect of the earth-flattening 
transformation on the ray path. 
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The spherical velocities are close to values given by 

Johnson (1967) for the upper mantle transition at 400 km. 

The transition in model X is assumed to begin at R = 6200 kn, 
i.e. below the low velocity zone. 


Model, WW. .of; d.. 6, km transition, at. 20 .km.im Fig. 4.1: 


Sy ey 6.63 km/s 
Sen 8.05 km/s 
og: = .2 

z = 20 km/s 247 2+ 0 km 

re) s r 

P = constant 

(4.50) 
Monel. of. 2 kmt transition at 20 km in Fig. 1.3: 

hae 6.63 km/s 
A a 8.05 km/s 
ge} Se * 025 

z = 20 km/s PODS Pal Oe ti 

re) s r 

9 = constant 


en 


Velocities used in models IV and V are characteristic for the 
Earth's crust and the models differ only in the thickness of 


the transition. 
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The density was kept constant in all the models 
which is not a realistic assumption. Basically the change 
of density must be very similar to change of velocity. However, 
the assumption about density being constant is not a serious 
limitation. The density variation introduces a small 
asymmetry into the model due to stretching of the vertical 
coordinate. This is a result of transformation (2.14) and 


provided the Birch's relation for density change is valid 


py Vasey + 1613 


(Wang, 1972), the effect of density variation on model X is 
shown in Fig. 4.18 and in Table 4.1. For models IV and V 
the distorted models are displayed in Tables 4.2 and 4.3, 
respectively. We see that the difference between the 
original and the distorted models is very small. 

For computation of synthetic seismograms the response 
integral (4.18) must be evaluated in wide range of frequencies. 
In a previous section the asymptotic pressure response was 
found in limited range of frequencies outside the limits of 
velidity of the aformentioned results, the numerical integra- 
tion must be done. This involves the calculation of the 
gamma function I'(z) and Hankel function 16) (2) for complex 
arguments. Programs for these were written with algorithms 


based on formulae valid in different regions of complex z plane 
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VELOCITY ( km/sec) 
9 9.5 LJ mS 


The effect of density variation on the 
velocity profile in Model X (full line). 
The triangles indicate the velocities in 
the model that includes density variation. 
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INFLUENCE OF DENSITY VARIATION ON 
VELOCITY PROFILE IN MODEL X 


Velocity Original Coordinate Transformed Coordinate 
v [km/s] z-z [km] c [km] 
8.41 -180 -169.6 
8.42 -160 ~15153 
8.43 -140 —132 08 
8.45 -120 —-124 2 
8.48 =100 = 9556 
8.54 - 80 = 26 99 
8.63 - 60 — §8 at 
3.28 - 40 -~ 39.1 
8.98 - 20 ~ 09 56 
R24 0 0 
B52 20 20.3 
2... 09 40 41.0 
10.00 60 @2Z 4 
LG. 03 80 S302 
EO) .20 100 105.4 
Le 120 Lied ee 
ae 9 AS fee) 140 149.3 
Le or 160 LlJviss 
10.38 180 193.43 


Table 4.1 
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INFLUENCE OF DENSITY VARIATION ON 
VELOCITY PROFILE IN MODEL IV 


Velocity Original Coordinate Transformed Coordinate 

v [km/s] z-z [km] co [km] 
6.63 -20. -18.69 
6.63 ot Se ee -14.02 
40:5 ~10. Ae < 
6.63 ae wont. 4 
6.64 =a = 0.95 
aa! - 0.5 - 0,49 
fae a 0. a 
(Ee Bi Geo the 
8.04 1 1.06 
8.05 Me 542 
8.05 10. 10.86 
8.05 5 LG. 32 
8.05 20% va Nae he 


Table 4/2 
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INFLUENCE OF DENSITY VARIATION ON 
VELOCITY PROFILE IN MODEL V 


Velocity Original Coordinate Transformed Coordinate 

v [km/s] z-z [km] C [km] 
6.63 —20. ~1s.0/ 
6..63 -15. -14.01 
6.63 -10. =~ 9.34 
6.63 Shy - 4.67 
6.63 abd. og — oS 
6.63 - “tou47 
6.65 - = 4,096 
6.76 -~— SOS - .049 
7.24 +U 4S 
fim 83 .05 5 Ode 
8.02 ook 2106 
§.05 vo ~54 
8.05 i. 25086 
8.05 Das 5 . Hs 
8.05 10. LO368 
8.05 Aids LG 
8.05 20 4 yA Ds he 


Table 4.3 
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(see Appendix D). It is easy to achieve accuracy comparable 
to that of other functions incorporated in the Fortran language. 
Thus the main concern is about the efficiency and 
accuracy of the integration procedure. The analytic solutions 
of Section 4.3 showed that the integration contour can be 
moved to pass through the saddle points via the steepest 
descent path. The approximations are only good at high 
frequencies and they deteriorate as frequency decreases. To 
find a suitable contour in the complex kK plane, the integrand 
in (4.18) must be studied. Its analytic behaviour is studied 
using contour maps of logarithm of its modulus and phase 
(Fig. 4.19). The grid depends on frequency, widths of saddle 
points and steepness of the descent paths. For example: in 
case of the triplication, different grids are used in regions 
near to caustics at 9. and 0. (Fig. <4. 29sand 4.91)% ithe 
signals at caustic 0, are weaker than signals at the caustic 
Ou which results narrower saddle points in the first case. 
As low frequency end is approached the saddle point 
corresponding due to the rays of branch AB disappears, this 
is due to the source being placed in inhomogeneous medium: 
the high frequency waves at shallow angles encounter the 
small change in velocity gradient around the source and con- 
tribute to signals forming branch AB. The long period waves 
are not effected and, therefore they contribute very little 


to amplitudes of branch AB. 
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When the steepest descent path is found on the 
contour map it is approximated by straight line segments 
(Fig. 4.19). On the contour map of the phase in Fig. 4.20, these 
approximately follow the constant phase line. Thus 
the numerical integration can be performed accurately using 
one of the simplest numerical methods - the Simpson's method. 
If the contour were not along a line of constant phase, the 
oscillatory character of the integrand would make the 
integration extremely difficult. According to the Simpson's 


method (Abramowitz and Stegun, 1965; p. 886) 


S541 
| 
f£(w,K,r)dk = BAK fe THOT “tHe + MK. Ft a Solel Paes ro 
K, 
J 
Cree} 
where f° fi ae f. are values of integrand at 


successive steps of interval Ta - Rigorously, the 


Lk hy 
j+1 
step size should be evaluated from the error term in the 
. 1 . . ° (iv) 
Simpson's method. This is proportional to f (K) because the 
Simpson's rule is correct to the second order, and it is rather 
tedious to evaluate the step size from it. Our integrand is 

a smoothly varying function along the path of integration and, 


therefore, the step size is simply estimated from the first 


derivative 
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AK 


| E(k)dk = Sk {£(-AK) + 4£(0) + £ (Ak) }-2— ft) EAR 
- AK 


ae 


where - AK < ey < Ak . If we assume that within (-Ak,Ak), 


the integrand is an exponential function: 


fie) = a * ALhe 


the fractional error 


Pee Ate hee 
Oo Go a See eS ea. 54 best 

180f(E,) Ak 180 
Thus [AK|.= 3.7 Eo Es anh (4.54) 

"Cr 
For a fractional error € = 10° we can take 
a, 7 #kKD — 
| Ak | ~ ITOET CR) Case 


where factor 3.7 was ignored for safety. Arg(Ak) is chosen 

in the direction of the contour segment. This method has been 
used by Phinney and Cathles (1969) and Chapman and Phinney 
(1972). The first derivative in the denominator can be 
estimated by taking the gradient at the endpoints of the 


segment. The step size of two neighbouring segments is 
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evaluated and finally their minimum is chosen for integration: 


[Ak (r)| 2 min{|Ak(«,,r)|, [AK(K, yor) |} (4.56) 
where 
Tun, .8.r)oK 
[AK (kK, 4r) | ” OTE CE, FORE) ~ CTC BE 
Then 
es S441 © 
j-1 
P(w,r,z_) = f rik,@,.hcr + 2 f f(K,w,r)dk + f tik oO, cak 
j=1 
= J K, Ky 


(4.57) 


Ki and K, are chosen so that the integrals to infinity (which 
are deep in valleys) are small. Again assuming exponential 


decay their contribution can be estimated as 


K 


; 7 £7 («) ,wsr) Cae E> 
{ r(K ,@, 4) ak? + j f£(K,wW,r)dK = —————_. - — 


J 


(4.58) 


The accuracy of the algorithm can be tested integrating around 


closed curves with or without singularities enclosed. All 


these trials proved to be within the required accuracy € = i 
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But, even if the integration is performed by 
Simpson's rule, the amount of the computing time needed 
would be enormous, if further savings were not made. Many 
ccemputations can be saved when one contour is used to 
evaluate the response at several ranges r . Only the Hankel 
function aS”) (cr) is dependent on r and only these values 
have to be recalculated for every r. The rest of the 
integrand is stored and reused every time the response is 
computed for new r . The contour of integration is no longer 
the approximation of the steepest descent path, but if chosen 
with care, a useful range of integrals can be evaluated with 
no significant loss in accuracy. Usually, the contour is 
chosen to be a steepest descent path for r which is in the 
middle of the interval of r's where the same contour is used. 
This is very efficient at low frequencies and in the neighbour- 
hood of the caustics where the saddle points change their 
position slowly with changing r . In the illuminated region, 
the contour must be changed more often as the saddle points 
separate more and more rapidly, until a separate contribution 
from each has to be computed. The phase difference between 
the two rays is growing and interference becomes important 
and causes oscillation of the total amplitude (Fig. 4.22). 
The contour integral cannot be evaluated along the real axis 


between the saddle points. It must be taken into the valley 
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above or below the real axis which is common to both saddle 
points (Fig. 4.23). The Figure 4.22 showing the results of 
the contour integration is similar to those shown and 
discussed in the former section. In the illuminated region 
the reflected wave of the reverse branch BC is stronger and 
the amplitude oscillates about its geometrical amplitude 
G4>30) . 

For the synthesis of pulses in the time domain the 
method described by Chapman and Phinney (1972) was used. In 
order to synthesize theoretical seismograms, spectral values 
over a wide range of frequencies are needed. To limit the 
Frequency range we have simulated a critically damped seismo- 
meter response with transfer function 

3 


Ww 


oS) bo eee Se (4.59) 
(9 -w*~ iw) (w2-w*-iw6) 


where Me and W, are natural frequencies of the seismometer and 
galvanometer. For critical damping, 46 = 2W and A = ce ; 
Although in period range used here Ue = .45 sec, Pe = 1.5 sec) 
seismometers do not use galvanometers, the transfer function 
(4.59) still represents a band limited instrument. 

The transfer function is convolved with the source 
of unit step function i.e. P6t) = H(t) . Their combined spectrum 


and the resultant seismogram are shown in Figures 4.24 and 4.25. 
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The combined spectrum of the critically 
damped seismometer and the unit.step 
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The response of the critically damped 
seismometer to the unit step function 
source. Marks on the time scale are 
every second. 





’ as< 
- - 7 
wa © 
yt yy » - 
» = o 


Se tees Tee : 





LZ 


The spectral amplitudes (4.29) and (4.30) are independent 

of frequency and, therefore, the pulses on forward branches 

of the travel-time curve have the same shape as in Figure 4.25. 
The pulses of the reversed branch will be shifted by - 5 id.e. 
they are "allied" pulses to that in Figure 4,25. 

The inverse Fourier transform (2.6) is evaluated 
using finite Fourier sum. Efficient programs, which use the 
algorithm for fast Fourier transform developed by Cooley and 
Tukey (1965), exist in Fortran and we only supply the spectral 
amplitudes. The spectral amplitudes calculated by evaluating 
the response integral (4.18) must be interpolated. P(w,r,0_) 
is complex and its modulus and phase must be interpolated 
separately. Where P(w,r,Z_) is evaluated, the phase is a 
discontinuous function in the range -7T to 7 . The geometrical 


results can be used to remove the 27 ambiguity in the phase 





between various frequencies and ranges. From (4.21) and (4.23) 
we get 
of ak wae dT 
or dr 
(4.60) 
of 
and Free T 
where f(k) is the continuous phase function, The equalities 


in (4.60) are valid only at high frequencies. At finite 


frequencies they can be replaced by approximate equalities. 
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As such they are used in the algorithm that forms continuous 
phase function needed for interpolation. In the region before 
the caustic the arrival time is formed by solving the stationary 
phase condition (4.22) in complex K plane. The arrival time 

is the real part of = fs) - The same applies to regions before 
the critical point where the transition is thin. In this case, 
however, the times of arrival of sub-critical reflections from 
a sharp velocity increase can also be used, especially where 
the transition is very thin. 

The theoretical seismograms evaluated using the 
Pranewer functions (4.59) with step function as a source are 
shown in Figures 4.26 to 4.29. They are normalized with 
respect to the maximum amplitude at caustic or at the critical 
point. 

The seismograms in Figures 4.26 and 4.27 for pressure 
at the near and far ends of the triplication. The maximum 
amplitude at 10°30' and at 20° is the Airy phase due to inter- 
ference. The geometrical endpoints are at 10°20' and 21° and 
as we move away from the illuminated region the signals decay 
rather quickly. The first pulse in both cases corresponds to 
the forward branches and is in phase with the source pulse. 

The pulses of the reverse branch BC are the allied pulses, 
shifted by - 7 with respect to the first arrivals. The 


2 


difference in phase can be clearly distinguished at ranges 
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NEAR END OF THE TRIPLICATION 
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Figure 4.26 The synthetic seismograms for Model X at 
ranges near the caustic C. The marks on the 
time scale are every second. The geometrical 
arrival times are indicated by longer vertical 
marks. 
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FAR}ENO OF THE TRIPLICATION - 10.4 K4/SEC VEL. IS USED TO REDUCE TIME 
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Figure 4.27 Synthetic seismograms for Model X at 
ranges near the caustic B. (Time scale 
as in Figure 4.26.) 
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Figure 4.28 Synthetic seismograms for Model IV at 
ranges near the critical point. (Time 
scale as in Figure 4.26.) 
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Figure 4.29 Synthetic seismograms for Model V at 
ranges near the critical point. (Time 
scale as in Figure 4.26.) 
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12°20' and 13° which are well inside the illuminated region. 
In this example, the response function (4.18) was calculated 


by numerical integration for 12 frequencies ranging from 


W =.5 rad/sec to w = 60 rad/sec. From w = 60 rad/sec to 
w = 110 rad/sec the geometrical values (4.29) and (4.30) 
were used. 17 spectral values were interpolated using cubic 


spline in order to obtain 128 spectral values between 


Ihe 


-82 rad/sec and w, = 104.7 rad/sec (Nyquist frequency) 


W 
N 


1 
with step Aw = .82 rad/sec. These were used as input to the 
inverse Fourier transformation. 

The seismograms in Figures 4.28 and 4.29 are for 
pressure in the neighbourhood of the critical point when the 
transition is thin. For both transitions the amplitudes are 


normalised with respect to amplitude at 70 km. For .2 km 


transition, this coincides with the maximum amplitude while 


for 1.6 km transition, maximum occurs at 80 km. This agrees 
with the shift of the critical points, Xo = 60.5 km and 67.8 km 
for .2 and 1.6 km transitions, respectively. In the region 


before the critical point, the seismograms show that partial 
reflection grows considerably as the transition becomes 
thinner. For .2 km transition, the signal at 20 km is quite 
strong while for 1.6 km transition it would barely show. Also 
their shape for the thinner transition appears to be more like 


the sub-critical reflection at sharp discontinuity. The partial 
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reflections are in phase with the source pulse which agrees 
with the fact that at sub-critical reflections, no phase change 
occurs. The amplitude maximum is due to interference of the 
reflected wave and head wave. Further beyond the critical 
range the two waves start to separate into first arriving 

head wave and later reflected wave. The reflected wave 

beyond the critical point is totally reflected and, therefore, 
shifted by - 5 with respect to the source pulse. This is 
clearly observable at 150 km. 

For calculation of the seismograms for thin 
transitions, the geometrical values (4.29) and (4.30) could 
not be used and the numerical evaluation of the pressure 
spectrum was performed at 17 frequencies between W = .5 rad/sec 
and w = 110 rad/sec. Their interpolation was performed in 
the same manner as for the triplication. 

The results presented in this section demonstrate 
that the numerical methods can be used efficiently to compute 
the spectral response and the synthetic seismograms. This 
should be an encouragement for their wider use in seismology. 
In application to the interpretation of data, the results 
for the response in the frequency domain are probably of 
greatest importance as the different properties are then 
easily seen and can be interpreted. Although results pre- 


sented here were obtained for relatively simple models, some 
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of their basic features (behaviour near caustics and critical 
points) should be very similar to those which one would get 
even for more complex models. Furthermore these results can 
serve as calibration results for approximate methods. Their 
importance has increased considerably in recent years but 

the error estimates have always been very vague. These can 
only be evaluated rigorously if the approximations are 


compared to exact results. 
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CHAPTER 5 DIFFRACTION AT THE VELOCITY REVERSAL 


Wave propagation in media with a monotonic velocity 
increase has been studied in the last chapter. The medium 
response was found and evaluated both exactly (or rather 
as accurately as is numerically required) and asymptotically 
(analytically). Triplications and critical points are not the 
only ambiguous regions for interpretation of travel-time 
curves. Problems also arise when a discontinuity is present 
such as the shadow due to the core-mantle interface or low 
velocity zone at the top of the mantle. The shadow is caused 
by a sudden decrease in velocity that may be abrupt or con- 
tinuous. This information cannot be found easily from the 
travel-time curves. The boundary between the illuminated 
region and shadow zone is not sharp even for a first order 
discontinuity in velocity structure. The diffracted signals 
propagating along the interface bring considerable amount of 
energy beyond the shadow boundary. In case of a continuous 
velocity decrease the amplitude of the signals decrease with 
distance so rapidly that a shadow is observed too. The 
amplitude decay of these signals into the shadow is very 
different from the decay into the "true" shadow. Only after 
amplitude studies can the nature of the velocity structure, 


which causes the shadow, be revealed. 


BB ba | 










; THT T HOTTOAATTIG 
sith I naoljeassvesq svaW 
cal ba lbut ased #64 ee 
bs siave baa bauvold é¢sv 98n01 as 
te: a 
< 
: r vi 19ML ef a6 visiatwI38 | 
. \ 
) > sti«c i »(vils 
2 ; a&god ;7 uy ugid 
al 7 1] oels .- idowzgd 
r s+ 5 - be 4 
> ) 74 a eb woOODe ug af 
3 io qo3 o03 3h Sonos 
3 to v mI sees toeb a 
} goon Ten olal wtf 
: ; paps: | 
. i 1or -39viusz smiz-Ie 
j7sfea jon el snos wobsde bas go 
- : ” 
‘ 5 : ; iSout3ia vitsoisyv ca! yjiuatinose 
> ae 
leo javoms sidatrebteanos anitxrd sossizejat sd3 grols abies 10 


evountinoo a -¥Ytebavod wobsde 







tw sensto9b 
La rola 


alsagts ofa to Pate = 
5! . es ‘ate ’ ian 
edt? , ae. 


' 








oda b hd xs 19, 
ie» 








died 





sony 





122 


The diffracted signal caused by a first order 
discontinuity has been studied by many authors (Scholte, 1956; 
Phinney and Alexander, 1966; Teng and Richards, 1969; Chapman 
and Phinney, 1972) although the interpretation of data is 
still uncertain (Phinney and Alexander, 19693; Bolt, 1970). 
This probably reflects regional differences in the core-mantle 
interface but may be caused by a velocity reversal just above 
foe interface (BokFt?2i970). No thorough treatment of the 
velocity reversal seems to have appeared in the literature. 

In particular the non-geometrical nature of the shadow and 
the decay rate of the signal does not appear to be fully 
appreciated (Bolt, 1972). 

The present chapter studies the wave propagation in 
Epstein medium with velocity maximum (Epstein, 1930). The 
results are compared to those obtained earlier for a parabolic 


profile and sharp velocity decrease. 


3.1 The reflection coefficient for velocity reversal 


Some of the properties mentioned in this section 
were studied previously by Epstein (1930), Phinney (1970) and 
their summary may be found in Budden (1961, p. 380) or 
Brekhovskikh (1960, p. 185). Some are included for the sake 
of completeness, others are essential for discussion of the 
response integral which none of the aformentioned authors 


attempted. 
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For Epstein profile with velocity maximum as and 


velocities above and below the transition vi (Sie.5 5.01), the 


vertical wave number (3.17) becomes 
2 2 
Ae ig (q, ~ qd, sech oe Core 


The reflection and transmission coefficients are 





Pia BS RS Ase 
P(1+2ioq,) P(s-2i0q,-y) ['(g-2igq,+y) , 
R ooo ERE ENS a le ee ee 
Sa F(-2icq,) Gay) T Gary) 
(one? 
EG 2k aye) 
2iog, I (s-2ioq,-y) I (4-2ioqg, ty) . beta Ake aie 
T(0.,0,,«) = ——=—_, + ______+ _ (5.3) 
iy (1-2ioq,) 
where for po constant 
Z 
2 VE {hts *5 
EL a oped: Se 
Y= {7 40 5 (1 5) 3 = {F Ss cos Pgh (5.4) 
Vv Vv 
1 oO 
may tess 
and OF tsofhe critical angle O. = sin age Instead of the 
oO 
absolute measure of the thickness of the layer o , the 
effective thickness § = al * 20 is introduced. The reflection 
ue 


and transmission coefficients depend on the thickness of the 
layer O as well as the wavelength of the propagating wave and 


the effective thickness S includes both effects. 
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VELOCITY REVERSAL 
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Figure 5.1 The travel-time curves and the rays 
propagating in the medium with 
velocity reversal. 
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The absolute values of the reflection and 


transmission coefficient are (Budden, 1961; p. 381) 


L. 


2 
a eosceny +*2 
[R | Pre 20y * cosh (2% S cos Ph 323) 
cosh (27 S cos 8) - 1 hs 
[tT] = (5.6) 


Serre ony +.c0s8h..(27;S.cos a)? 


For thin layers or long wavelengths the effective thickness 


is small and the reflection is very weak. 


Sey ks and |R| + 0 


Pol em 


For thick layers or short wavelengths when the effective 
thickness is large the amount of reflection depends on the 


emergent angle from the source 61° (Brekhovskikh, 1960; p. 185) 
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At angles oy less than critical most energy is transmitted 
while at angles of emergence greater than O. most energy 
is reflected. The range of angles when both, the reflection 
and the transmission, occur depends on the effective thick- 
ness (Fig. 5.2). The boundary between total transmission 
and total reflection is blurred by partial reflection at 
sub-critical angles and tunnelling effects at angles beyond 
critical. Both partial reflection and transmission by 


tunnelling increase as the effective thickness decrease. 


-7™S(cos O,-cos OM 
[R| we ; oe. =< 6 


= a aS 
TS (cos o cos 01) 


|T| me Os eo 


At this point it is desirable to point out the advantage 
of the Epstein profile over the parabolic profile in 
Fig. 3.3. Ihe solutions of thé wave equation” for a para- 


bolic velocity structure 
Oo 
VY = vert oe ie (5,8) 


can be expressed in terms of parabolic cylinder functions 


(Phinney, 1970). Then the reflection coefficient is 


|r| = (5.9) 
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Figure 5.2 The dependence of the modulus of the 
reflection coefficient on the effec- 
tive thickness of the "sech'" layer 


for v, = v, = 7.83 km/s, v_ = 8.05 km/s. 
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Figure 5.3 Parabolic velocity profile; 


a= |Fal, 
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where a* is a complex conjugate of 


2,- Oe 2.> 
Se Ga fe) _ 2 Vv 
a ae ce - (<8) (sin? & - vy C3)) (5.10) 
fe) re) Vv 
fe) 
. +s 
and oP = l= | is so-called half-thickness of the 
oO 
parabolic layer. The reflection coefficient (5.9) was first 


deduced by Rydbeck (1943) for propagation of waves in the 


ionosphere. It is evident that it does not give correct 
oO 
answers for thin layers (5¢ >i0)\which is R 0 rather 
1 
than R %™ — . This discrepancy is caused by discontinuities 
2 


of the velocity gradient on the boundaries of the layer. 
Northover (1962) found a power series solution for thin 
parabolic layers which yields the correct reflection 
coefficient. His results were recently confirmed by Chesell 
(1971) who determined the reflection coefficient by exact 
numerical integration of the wave equation (Fig. 5.4). Northover 
(1962) criticizes the method of matching the solution above 

and below the transition zone in order to satisfy the boundary 
conditions. While partly justified in the case of the parabolic 
layer solved by Rydbeck (1943), his criticism does not apply 

to the Epstein models. The latter has no gradient dis- 
continuity and no additional conditions have to be satisfied. 
The correct values of the reflection coefficient (5.5) at 


different limiting cases prove that the method of matching 
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Figure 5.4 Variation of the modulus of the 
reflection coefficient with frequency 
for parabolic layers of different 
effective thicknessess uU = sh Pee 


f and A_ are the critical frequency 
and wavelength of the ionospheric 
layer (from Chessell, 1971). 
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the solutions may be used in Epstein's problem. It can only 
be used in a limited range of frequencies and thicknesses for 
the parabolic layer. Subject to the aformentioned constraints 
the reflection coefficient (5.9) gives useful results. 

In the complex kK plane the reflection coefficient 
(5.2) has three sets of poles due to the gamma functions and 


. W) U) Ww 
Pec ered points th. = f—~ 2k | = t-—— and ti€ = #— . 
1 Vi s ae = sf 


Of total number eight only two of the Riemann surfaces are 
important in our discussions. The branch cuts from s to 
ry, wteeisctnity, woe ly s oF fr, are very close to each other 
and the top (+) and the bottom (-) sheets of qd, are taken in 
agreement with the convention introduced in Section 4.1. 

The secular equations that determine the position 


of the poles are 


2ioq, = -n rie a Ek Bie Ske See es (a) 

2ioq, =n+%+#y¥ 1 ek | bls RAR es Tiare (b) 

2ioq, =n+%- y navd®,Le2Z, Bvene Ge) 
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The sets of poles are given by (5.1la) 
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are on the top Riemann sheet "beyond" the branch points 


Ww 
+—., (Fig. 5.5a). The signals arising from the contri- 


Vv 
1 
bution of these poles propagate in horizontal direction and 


decay exponentially with depth. Their phase velocity 


L 
-* 


is smaller than vi and they arrive to the receiver much 
later than signals we are interested in. 


The sets of poles given by (5.11b) are 


ho 
wer 


+ (n + 3 + y) 7/407] Ge ig 
Z 

At low frequencies when y is real the poles lie on the 

real axis "beyond" the branch points + ~- on the lower 


v 
1 
Riemann sheet (Fig. 5.6b). The zero frequency position is 


oO n+1 
2h 2G 





(5.14) 


As frequency increases they move towards the branch points. 


At frequency w* 
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2 5 
W 2 2 
* = es 1. 
+ ks a 5 + (n + 43) /40 ] C56.) 
aff 


and they leave the real axis. As frequency increases they 

move anticlockwise in the first and third quadrants of the lower 
Riemann sheet, their real part decreasing in absolute value 

and imaginary part increasing in absolute value. They reach 

the branch cut where |Re,k | v a and pass onto the top 

Riemann sheet. Then they continue to move towards the high 


frequency position where they slow down and "settle" at 


2 lL. 
2 1 2 
ak, = t + f(nts)? + 3/407 + a (nts) $5 (cos 6-5) 7] 
Vv 40 4S 
Oo 
pat, 


(Fig. 5.5a). The complex poles on the top Riemann sheet 
contribute to signals propagating with phase velocity between 
vi and Lae and are important in the evaluation of the response 
integral. Only the poles in the first quadrant contribute, 
however, as they are close to the steepest descent path. 
Further discussion of the behaviour of the reflection coeffi- 
cient in the neighbourhood of these poles is left to the 


next section. 


The two sets of poles given by (5.llc) are 
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and they lie on the lower Riemann sheet. Their zero frequency 
positions are 


+ K = + 


n 
an 7 4 26 (5.19) 


As frequency increases the poles move along the real axis 

away from the branch points. As w = w* they reach the 
position oKF and coalesce with the set of poles aKa Gs Te Kes 
They leave the real axis in the direction apposite to that 

of aKa and move clockwise in the second and fourth quadrants of 
the lower Riemann sheet, (Fig. 5.6b). As high frequency 


limit is approached, their movement slows down and they 


settle at positions 





2 ks 

aku 7 * Fg t+ Lente)” + S15 - ants) 45 (eos? 0 - 2’ 
2 Le 40 40 4s 

(5.20) 


(fie. 3.5b). - The. poles 32 » being on the bottom sheet, 


have no effect on the reflection coefficient and no signals 


arise from their contribution. 
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5.2 The response integral 


The rays most effected by the velocity reversal are 
those emerging at angles close to critical angle o. «) “Lee 
response integral (2.22) for the velocity reversal where 


density is constant is 


+o 
P(w,r,t_) = + P (ww) f R(G 655K) BS) Cer) i dk 


a ga Og 


where R(T 555K) is the reflection coefficient given by (5.2). 
At high frequencies the same method of evaluation can be used 
as that used previously for the monotonic velocity transition 
(see Section 4.3). The reflection coefficient can be approxi- 


mated by the WKBJ reflection coefficient (4.17) and its phase 


is 
mage 1 . 1 
404) (4-2i0q, +7) (s-y) 
o(K) = 20q, log ————————.— - iy log A tC Se ee” 
aoe tae eG CS. c280 


The point of stationary phase of the integrand kK* is given by 


the condition 
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The equation (5.23) must be solved numerically for complex kK 
_ This gives two saddle points very close to the real K axis 

W W , > : 
between K = — and k = — if the horizontal distance is 

fe) Vv ts vi 
greater than Xo 2 There are two arrivals at these distances 
(Fig. 5.1) which at high frequencies are represented by two 
aforementioned saddle points. The saddle point close to — 
oO 

represents the ray close to the critical ray and is most 
effected by the velocity reversal while the other represents 


the ray reflected by the large velocity gradient (Fig. 5.5a). 


The amplitude at the saddle point K* is proportional 


to 
dr ty: d*6 -% 
pa Pilbes ailid roe aha Crees) 
K* ak ik* 
where 
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At high frequencies the saddle point approximation gives the 


response 


B ,(w) 2 % ike + id(K*) 


BGO thet shoe beet te ecores eavrie €5',.37) 
2/27 





for the wave propagating along the reversal. The saddle point 


approximation for the reflected wave is 


Pp lk iKkk*r + j * 
P(w,r,t_) “8 21 —— aici Waite A (5.28) 
‘i suee {-r ° aay a 
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The geometrical arrival time is given by 


al 
il 
Ele 


Re{k*r + o(k*) } (5529) 


The second order saddle point approximation fails near the 
ehustdce€r(Picuree5.1)4 In its vicinity the same treatment 

as that described in detail in Section 4.3 can be used to 

obtain the high frequency approximation in terms of Airy 
function. The behaviour in the caustic region was investigated 
in Chapter 4 while in this chapter the attention is concentrated 
on the effect of the velocity maximum on the wave propagation. 
As the angle of emergence 9. > Ps approaches the critical angle 


1 


the ray travels further in the region of the reversal and its 
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amplitude decays. The saddle point which is the approximation 
of this ray at high frequencies approaches the zeroth pole 

in the string, Kn (5.13). The contribution from the integra- 
tion through the saddle point becomes equivalent to that of 
the residue of the zeroth pole (Fig. 5.7). The amplitude 


given by the residue of the zeroth order pole decays exponen- 


tially when r is large (Appendix C). 


* ~ B(w) 116 fa eA asi oe On Mocks uce 
ee. & ) 4 = (= ° 
r 40 q 
So Ker 
2/0 
Llo sty)? + 2 1 (5-7) (5 30) 
aa ifY ae ap re 
w? 2 2 
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The decay rate of the signal is given by the imaginary part 


of the zeroth pole. To first order in frequency 


2 1 
Vv Vv *, 
Im .K_ = (n+) mee PS Clive —t) pep ke 3 C5 391) 
etal 20 2 y 
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Thus the decay rate at high frequencies is independent of 
Frequency, The dashed line in Figure 5.8 represents the 
residue evaluation of the response integral (5.21). When 
the numerical integration was computed for this signal 
separately, the amplitude curves coincided with this line at 


frequencies w 2 8 for distances X = r 2 700 km. 
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Figure 5.7 The steepest descent path in the complex 
K plane for large horizontal distances. 
The crosses indicate the poles 9%? 
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The amplitude curves evaluated numerically for 
different frequencies (full lines). The geo- 
metrical amplitudes from the saddle point 


evaluation (chained lines) and the residue of 


the zeroth pole (broken line) are also indicated. 
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Chapman (1969) studied the decay rate of rays 
effected by the parabolic velocity reversal (5.8). Then the 


reflection coefficient has poles at 


2 ail oe 
eee ec es |—-| 
n Vv Vv 
Oo oO 


(n+%5) S42) 


and the decay of the low order poles 


vv L 
-r|—"| (n+) 
u% @ “ 
also independent of frequency. Both profiles have the same 


curvature at the maximum if 


1 ve 
AS 5 ey. = 1) (5.34) 
v 4o Vv 
fo) 1 


and the decay of the amplitudes with growing distance is 
identical then. This result could be anticipated as these 
rays are effected mainly by the velocity structure at the 


maximum. Chapman (personal communication) showed, that the 
-l 
dr “ 
geometrical amplitude given by lr Wal 


reduces to the same 
exponential decay as that given by the zeroth order pole. 

It can be said that due to the rapid decay of the amplitudes 

a "shadow" is observed. Although the signals are present their 
amplitude is so small that they cannot be detected. Thus the 


"shadow" caused by the velocity reversal is very different 


from the real shadow caused by an interface when no rays 
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propagate beyond the shadow boundary. 

The signals near a real shadow caused by an inter- 
face are diffracted signals grazing the interface. These 
diffracted signals propagate along the interface and emerge 
at distances beyond the shadow boundary given by the grazing 
ray and, therefore, blur the shadow boundary. Their contribu- 
tion is obtained from the residues at the poles of the coeffi- 
cient of reflection from the interface. Duwalo and Jacobs 
(1959) found the approximate positions of the low order poles 
for a model of a fluid sphere in a homogeneous medium. Their 


result (in spherical coordinates R, 9, o>) 


2/3 it 
ot he 37 3 
= 1. ———_ 
‘ee KOR. or ah ap. kc) (Tk R ) Ew. 
ie 
: Lo ; 
gives a decay constant proportional to w - (v is the wave- 


number in O direction i.e. equivalent to * in our case; R. 


is the radius of the sphere and ky is the wavenumber 


ee 
a(R) 
at the surface of the sphere.) The position of poles is 
slightly modified if a layered or inhomogeneous medium is 
assumed (Phinney and Alexander, 1966, 1969, respectively) 
but the asymptotic frequency dependence remains unchanged. 
The constant c depends on the boundary conditions and for 
the models just mentioned |c] < 1. It is, however, rather 


difficult to find c analytically. Chapman (1969) also gives 


useful approximate expressions for positions of poles when 
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the parabolic velocity reversal lies near an interface. Then 
the decay rate depends on the relative position of the inter- 
face and the reversal and on the frequency range considered. 

The diffracted signals propagating along the inter- 
face or velocity reversal are mathematically given by contri- 
bution from poles of the reflection coefficient. The dis- 
continuity on the travel-time curve due to the shadow can 
only be interpreted correctly when the amplitudes are analysed. 
The frequency dependent amplitude decay is basically different 
from the frequency independent decay due to the reversal. 

At short distances the evaluation of the integral 
by method of residues is not valid (Appendix C) and the high 
frequency approximation is given by (5.27) indicated by the 
chained line in Figure 5.8. This can be called the "illumina- 
ted" region in analogy to the case of an interface shadow 
for which the residue evaluation is not valid either before 
the shadow boundary. The region when the saddle point 
evaluation (5.27) is not a good approximation due to the 
poles’ influence, and the high frequency approximation is 
given by residues (5.30), can be considered a "shadow". As 
frequency decreases the poles' positions become frequency 
dependent and both approximations become invalid and the 
response integral must be evaluated numerically. Its values 


are shown in Figure 5.8 in dependence on range r = X . The 
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oscillatory character of the total amplitude for w = 8 is 
given by the interference of the two waves arriving at the 
receiver almost together near the cusp C (Fig. 5.1). The 
amplitude oscillates around the geometrical values (5.28) 

of the wave reflected by the velocity gradient. These waves 
form the reversed branch of the travel-time curve in Figure 
5.1 and are much stronger than the waves propagating along 

the reversal whose geometrical amplitude at large ranges is 
given by the pole's contribution (5.30). These come to the 
receiver earlier than the reflected waves and form the forward 
branch CD. All the amplitudes for frequencies larger than 

W = 8 oscillate around the geometrical value of the reflected 
wave. For frequencies w < 3 the total amplitude is smaller 
than the geometrical value and decreases rather rapidly with 
frequency. The transmission of the energy through the velocity 
barrier increases at low frequencies and the amount of reflec- 
tion decreases, 

The results of this’ ‘chapter inditate that “diffrac- 
tion" at velocity reversals depends entirely on the velocity 
structure near the maximum. The "shadow" is caused by large 
geometrical spreading along the ray path. Using the geometrical 
ray theory, Chapman (personal communication) has shown that 
its growth is exponential. This is confirmed by the results 


of the full wave theory presented in this section. The 
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equivalence of the parabolic and es 


profiles is a very 
encouraging conclusion as it gives us a tool for investigation 
of shadow caused by more complex structures. Provided the 


a 


profile has the same curvature as the part of the 
complicated velocity structure that causes the shadow, the 


amplitude decay should be equivalent. 


5.3 Synthetic seismograms 


The parameters of the Model II which was used to 
compute the amplitude curves in Figure 5.8 are: 


Vv. oF a 7.63 km/s 


v = 8.05 km/s 
re) 


o = 10 (5.34) 
2 itm 35 km a = 0 km 

fe) s r 

O = constant 


The velocities were chosen close to velocities at the top of 
the low velocity zone in the Earth's mantle. 

Theoretical seismograms in Figures 5.9 and 5.10 were 
computed using the method described in Section 4.4. First the 
response integral (5.21) was calculated by numerical contour 
integration for 8 frequencies ranging from w = .5 rad/sec to 
w = 30 rad/sec. From w = 60 rad/sec to w = 110 rad/sec the 


high frequency approximations were used. For the reflected 
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10.KM/SEC VEL. IS USED TO REDUCE TIME 
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The synthetic seismograms for the Model II normalized 
with respect to the Airy phase amplitude at X = 500 km. 
Marks on-the time scale are every second and the 
geometrical arrival times are indicated by longer 
vertical marks. 
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SHADOW FROM VELOCISY REVERSAL - 10.KM/SEC VEL. IS USED TO REDUCE TIME 
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The synthetic seismograms for the Model II 
normalized with respect to the amplitude 
at X = 700 km. (Time scale as in Figure 
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wave (5.28) was calculated for all distances. The amplitude 
of the "diffracted" wave (or wave propagating along the 
reversal) was approximated by (5.227) when r = xX < 650 km 
ana by (5.30) when r.= x 3% 650 km. 17 spectral values were 
interpolated using cubic spline to obtain the complete 
Spectrum (Section 4.4). The unit step function source 
P_(t) = H(t) and the transfer function of a seismometer 
given by (4.59) was assumed. Their convolution gives the 
spectrum in Figure 4.23 and the response of the seismometer 
to P(t) is in Figure_4.24. 

The synthetic seismograms in Figures 5.9 and 5.10 
show clearly the rapid decay of the first arriving wave 
that propagates along the maximum. The seismograms in 
Figure 5.9 are normalized with respect to the maximum ampli- 
tude at X = 500 km given by the Airy phase due to caustic 
C in Figure 5.1. The marks on the time scale are every second 
and the geometrical arrivals are marked by longer vertical 
marks. The first arrival at X = 900 km is negligable when 
normalized with respect to the amplitude at the catatice For 
greater distances the first arriving signal cannot be seen. 
The same seismograms but normalized with respect to the 
amplitude at x = 700 km are shown in Figure 5.10. The second 
arrivals are the waves reflected at the velocity gradient and 
are much stronger. At large distances only these can be 


observed. 
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CHAPTER 6 CONVERGENCE OF THE RAY EXPANSION 


Theoretical seismograms presented in this thesis 
were computed from exact Epstein's solutions (3.20) of the 
wave equation (Epstein, 1930). One question arises immediately: 
what would the seismograms look like if the Epstein profiles 
were replaced by a stratified medium consisting of homogeneous 
layers. The last has been studied extensively by many authors 
(Thomson, 1950; Haskell, 1953; Brekhovskikh, 1960; Gilbert 
and Backus, 1966; Muller, 1968, 1969; Hron, F. and Kanasewich, 1971) 
as it represents the simplest approach to studies of wave 
propagation within the Earth. Cisternas et al. (1973) 
showed that the exact solution of this problem, represented 
by Thomson-Haskell matrices can be expressed as infinite series 
whose terms can be physically interpreted as rays. Each ray 
is represented by a product of reflection and transmission 
coefficients at the interfaces on its path. These are inde- 
pendent of frequency and Cagniard-de Hoop can be applied to 
obtain the solutions in the time domain for each ray (Cagniard, 
£902;°02 Hoop, 1960) “ir the contribution from all trays in 
the "complete" ray expansion given by Cisternas et al. (1973) 
were evaluated and added to give the exact response of the 
medium, enormous amount of computing time would be required 
even for a small number of layers. Thus only the most important 


rays, i.e. those of maximum amplitude, are usually included 
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and the rest of the series is neglected. The use of the 
"partial" ray expansion is only justified if the ray series 
is convergent. The convergence, while assumed, has never 
been proved analytically, and the reason probably is that 
it is rather difficult. Some of the problems encountered 
when attempting to prove the convergence of the ray series 
as well.as a proof in one special case are discussed in 
Section 6.3. 

The wave equation can be integrated numerically 
for the continuous velocity profile and for the layered 
medium. The solution to the first problem is equivalent 
to the analytic solution obtained by Epstein while the other 
is equivalent to the complete ray expansion. The errors 
introduced by the approximation of the continuous model by 
layered medium are studied in Section 6.2. We also study 
the partial ray expansion that includes only the once reflected 
rays and investigate the error arising from neglecting rays 


with more reflections. 


6.1 Fundamental solution 


Equations of motion (2.4) and the constitutive 


relation for pressure (2.3) can be written in the form 


2 
u 0 - ices SA 2 ae 


2 
aS = & fz) (6.1) 
P pu” 0 ¥ 
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after the transformations with respect to time t andr. 


The equations (6.1) can be written in the matrix form 


dV (z) 
merece BLS F, Ch. 2) 


and the solution of such an equation is given by 


W(z)\ = P(z,2.) Vn) (6.3) 


where P(z,2.) is so-called propagator, introduced to 
seismology by Gilbert and Backus (1966). It is a fundamental 
matrix such that P(z522,) =e (identity matrix) and it can 


be found from 


PCase J Ease (z 0) (6.4) 


where F(z) is any fundamental matrix. A fundamental matrix 
is a non-singular solution of (6.2). Gilbert and Backus 


(1966) stated that formally 


z 
P(z,2,) = exp | M(E) dé (6.5) 
Z 


which is applicable particularly if M(&) is independent of z. 
Thus, if we assume a medium of homogeneous layers (Fig. 6.1) 
and choose an intermediate point ey in i-th layer 


M(z) = M(E,) 254 Pe Zs (6.64) 
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Figure 6.1 The stratified medium of homogeneous 
layers. 
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The propagator (6.5) is 


P(z,52;_4) = € (6.7) 


where M. = MCE.) and dy =eZy <2 


7 i i-1° We denote N Asmatrizx 


of eigenvectors of MN ewine. 


MsN-=N'B (6.8) 


Co) 


where B is a diagonal matrix of eigenvalues (we use the 


notaticn introduced by Chapman (1973)). The columns of N 
Co) 


are plane wave solutions and the elements of B 


pS a yt. N (6.9) 
; Z fo). 2 (9) 
are their wavenumbers. We denote Ny = NCE5) and By B (€,) 
where i corresponds to the i-th homogeneous layer. The 
propagation in the i-th layer is given by 
TE ae eer s 
Meets eye wae Oo Jee ed eo eo Be. ok ey N, 
pes Boi Ro = =i i =ii =i =i 
(6.10) 
since HS cy San Cee Pane 
hey = =] SL 
Thus the solution at the i-th layer can be written 
ar -1 
V(z,) = N, e N. ° V(z,_4) (6,23) 
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and the physical interpretation is evident: ea resolves the 


=i 
displacement and stress at z = 254 into plane wave components 
2a, 
© multiplies each component by its phase across the 


i-th layer and N; recombines the plane waves into displace- 
ment and stress at z = Ze The total propagator across the 


layered medium was given by Gilbert and Backus (1966) 


K 
P(Zy52,) = Tl P(z,525_4) (6o12) 
i=l 
Therefore it contains products 
eee (6.13) 
=i =itl =i ; 


which represent the transfer function of the i-th interface. 
The vector of amplitudes of individual waves at 


the (i-1l1)-th interface is 


) (6.14) 


while the vector of complete solutions in terms of plane 


waves including phase function in the i-th layer is 


) (6.15) 


The reflection-transmission problem at the i-th interface is 


a solution of 
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WG Es oe) (= Vz, + 0) (6.16) 


a" By XA) * Naa 2 


There will be as many possible reflection and transmission 
problems as there are columns of N (see Figure 3.5 in Chapter 
3 for acoustic case). Thus, the reflection-transmission 
problem can be solved using matrices at each interface, 
which, of course, is nothing new because it represents the 
Thomson-Haskell method, developed twenty years ago (Thomson, 
1950; Haskell, 1953). Thus the Thomson-Haskell matrix method 


is equivalent to finding the propagator P(Z, 92.) 


pg Big Bq 
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1g = I Be 
and V(z,) = P(Zy52,) ba Oe (6.0155 
Cisternas ét.al. (1973) showed that if.a vector xX 
is formed of all vectors x, Lii@ tls 2b 4d ee CERO TEL lLeciion-= 


transmission problem in a layered medium can be written in the 


following form 


X=QX+X (6.19) 
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where Xo is the source term. This can be expanded into 


infinite series 
co 
Some eek (6.20) 


and Cisternas et al. (1973) verified that terms in this 
power series represent physical rays. This is an important 
conclusion because it gives a new interpretation to the 
Thomson-Haskell matrix solution which is a complete solution 
Sey (652). 

Hron, F. (1972) showed how to generate unconverted 
rays in such layered medium and the formulae he obtained 
determine all rays possible in the layered medium. He 
introduced to seismology so-called "kinematic" and "dynamic 
analogues" in order to describe groups of rays with identical 
properties. The waves which travel in the layered medium 
along different paths but with identical travel-time curves 
(Fig. 6.2) are kinematically equivalent and thus called 
"kinematic analogues". Their total number Ny depends on the 


half-number of segments in each layer n, (Hrons F.; 12971) 


J-l n,+n, -1 
N.(n,; re n,) = eee ip itl C6521) 
i=1 “441 
sobishaha n,+n -1 (a. + ft = LL) 
Peg Ol: 2 eMart i+1 ; 
n ree +B. 
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are combination numbers. The group of kinematic analogues 
consists of several subgroups of waves which have identical 
amplitudes. These are called groups of "dynamic analogues" 
(Fig. 6.3) and the number of waves within each group depends 


on the number of reflections m at each interface: 


J-1 n, M4471 
N,(n,, ree TTS Mis cee my_4) = iT Ja a -] (6.22) 
i=1 Pe i 
where Vs = ms is the number of transmissions through the 


d-th interface (Hron,>F., 1971). The amplitude of the whole 
group of dynamic analogues called sometimes "total dynamic 
effect of the group" is evaluated as a product of the 


amplitude of one dynamic analogue and the number N This 


4 
saves a great deal of computations where the theoretical 
seismograms are calculated. 

Based on results by Hron, F. (1972), Chapman (personal 
communication) showed that the complete unconverted ray 
expansion in the medium consisting of K homogeneous layers 
overlaying a homogeneous half-space (Figure 6.1) is given by 
the expression in Figure 6.4. Although this expansion is 
basically equivalent to that of Cisternas et al. its terms 
can be interpreted more easily than the complicated products 


of matrices in the series (6.20). The time transformation 


given by the Cagniard-de Hoop method is (Helmberger, 1968) 
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J 
T = pr + : 2n ny d CG.ico) 


and the term 


n,-1 

n -1 

J-1 
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my_,7max(0,n,_,-n,) (6.24) 


represents the total amplitude of all rays with identical 
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al 4°22 
reflection from the top and bottom of the i-th interface, 
are the coefficients of trans- 


respectively. and 
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mission through the i-th interface in positive and negative 
z direction, respectively (Figure 6.1). 

In practise only the "partial" ray expansion is 
used, which means that for given number of layers K only a 
finite number of kinematic groups is considered, and within 
these only the dynamic groups with significant dynamic effect 
are taken into account (Hron, F. and Kanasewich, 1971; 


Muller, 1970; Gilbert and Helmberger, 1972). The partial 
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ray expansion may be viewed as an approximate solution of 
(6.2) that does not satisfy exactly the boundary conditions 
at each interface. 

The equation (6.2) can be solved numerically for 
both models, the Epstein monotonic transition and a stratified 
medium of homogeneous layers. The first yields the exact 
solutions (3.20) and the other will give the complete ray 
expansion in Figure 6.4 which should converge to the first 
one as the layer thickness decreases. The partial ray 
expansion can be obtained adding rays with certain number 
of reflections at each interface. The behaviour of the 
solutions with respect to the layer thicknesses is studied 


im Che wanext section. 


6.2 The numerical solution of the wave equation 


The components of the vector V(z) and the wave- 
numbers K and Q are complex functions and, therefore, 
equation (6.2) represents four equations. The equation (6.2) 
can be numerically integrated step by step as Q = Q(z) is 
a knowr function. Many numerical methods exist for integra- 
tion of systems such as (6.1). A very efficient method by 
Bulirsch and Stoer (1966) was used to obtain the solution 
¥V¢z). For this method, the initial values must be known and 


the accuracy of the integration depends on their accuracy. 
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In our problem it is not difficult to determine the initial 
values as in regions far above or far below the anomalous 
region centered around Zaz the plane wave solutions are 
valid (see Chapters 2 and 3). 


For the monotonic velocity increase given by 


— 45 
v(t) = (Sp 5 e'/%) Cr + eF/%) (6.25) 
where © = (z-z_) Laie? Sees. .) £CCOLGLNE tO a. )” 


the vertical wavenumber is 


Q(t) = ae eae - kK J (6.26) 


Beever ene TUrNning point. ot a ray, ¢ two travelling solutions 


ie 
should be present while below the turning point only the 


exponentially decaying solution is valid. From this, the 


initial values for the integration can be set as 
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K remains constant and represents the parameter of the ray 
whose turning point is at Co (4.24). After several tests of 
the program it was found that the initial values should be 
taken at the depth Cp which is approximately four wave- 
lengths below the transition. This seems to be an optimum 
depth with respect to the accuracy of the results and the 
computing time required to calculate them. 

le is evident fromythe initial values (6.23) and 
the matrix M in (6.1) that the numerical solution WUC of 
(6.2) remains real for all t . Thus above the turning point 
it represents a standing wave rather than travelling waves 
we would like to study. Analogously to (6.14) these solutions 


can be obtained from the matrix 


We 


X(c) = N oN CE) (6.28) 


The matrix of eigenvectors N is given by 


iQ(t) SAH tz) 
2 A! 
OW OW 
ee Ie (6.29) 
} 1 


and the vector X(S) is formed by the solutions representing 
the waves travelling in the positive (X,) and negative (X_) 


C-direction: 
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xX, \ 
X(t) = (6.30) 
a 
; -l 
The inverse matrix N is 
2 
gs oe 
2Q(7) 
-1 
N"(S) = (6.31) 
2 
10W 
2Q(S) 7 
and from (6.1) we obtain 
wp ad See 
= aP TO ee 
if 3 Bessie 98 
wey ipw 
a 6P + 20 ur 


where i and P are components of V(t) determined by the 
numerical integration. The reflection coefficient being the 


ratio of the upgoing and the downgoing waves is 
R =e (62453 


where ‘ 8 = arctan 





The numerical solutions of (6.1) for Model IV (4.50) 
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denoted Xo and X . These are compared with numerical 
solutions of (6.1) for a stratified medium of homogeneous 
layers which are denoted ee and _X . ‘The last represent 
the complete ray expansion in Figure 6.4 and should be better 
approximations of x and X as number of layers increases. 
The layered medium is obtained from the continuous velocity 
“Model IV (4.50) by dividing it into equal steps in velocity 
6v . The modulus of x, is compared at the point ce and at 
the turning point or to the modulus of the complete ray 
expansion Frese The point Gs lies 20 km above the center 

of the transition i.e. at the same level as the source and 
receiver for synthetic seismograms in Chapter 4. The percen- 


tage error 


a te 


|X, 

of Model IV, w = 100 rad/sec and’) ~ .5 km is givenGn 

wante, 6... It is evident that at this frequency the results 
for both models become equivalent only if the number of layers 
is very large. The percentage errors at both points are less 
than 2% if the number of layers is greater than 400. That is 
such a great number of layers that if the Cagniard-de Hoop 
method were used to evaluate the synthetic seismograms it 


would consume enormous amount of computing time. Thus we can 
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THE PERCENTAGE ERROR OF THE COMPLETE RAY EXPANSION 
FOR MODEL IV AND w = 100 rad/sec 





Number of Layers Beror -2t C-7720.km Error 2t or 
K E[%] e[%] 
3 99 98 99 9% 
4 98.96 oo whe 
5 ve ag ppg 
7 84512 83.8] 
10 68.44 bd wd S 
1 fhe 30205 14 a Ba 
20 20,03 i ho Sa? 
30 26.18 261.99 
40 19.84 20.41 
50 16.05 16234 
60 £3390 L348 
70 11.66 Wer? 8 
80 10.24 bat Pay 
90 9.108 S.09) 
100 8.186 8.144 
120 6.766 6.83 
140 I Pete: See tls 
160 4.946 4.940 
180 4.336 4.347 
200 3.853 3.963 
250 3.005 SP ALN Bs 
300 2.465 Fe}: fi 
350 2.095 are Ug 
400 1.828 i 2870 
450 LP ae oS es 
500 1.465 1.479 
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again appreciate the value of the exact solution (3.20) from 
which we were able to compute the exact synthetic seismograms 
without running into the aforementioned difficulty. Muller 


(1970) computed theoretical seismograms for a linear transi- 


tion 1. km thick with velocities v 6.4 km/s, v 8.2 km/s 


us 2 


an. 3 beat coms - He found 


i 


and densities she 5 0 bid canis and P, 
that the seismograms at epicentral distance of 150 km for 
the linear transition are equivalent to seismograms computed 


for a layered medium with 240 layers. The transition in our 
Vv 


Model IV (4.50) is thicker and continuous and the ratio a = .,82 
2 
is greater. Thus the first guess on the basis of Muller's 


conclusion is that the number of layers necessary for a good 
approximation of Model IV will be even greater. These 
qualitative estimates are rigorously justified by our 
quantitative results presented in Table 6.1. 

We have already mentioned that the synthetic seismo- 
grams if evaluated by Cagniard-de Hoop method require a great 
deal of computing time. That is why only the most important 
rays i.e. rays whose amplitude are the largest, are considered, 
when the computation of the synthetic seismogram is performed. 
The amplitude is dependent on number of reflections and 
transmissions the ray suffers on its path from the source to 
the receiver. We will study the convergence of the ray 


expansion analytically in the next section. In this section, 
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however, we would like to demonstrate some of the weaknesses 
that the partial ray expansion can have. The accuracy of 
the partial ray expansion depends on the correct choice of 
rays that contribute most to the amplitude. Such choice is 
often biased by the personal judgement of the authors - usually 
the rays that suffer only one reflection are considered to 
be of greatest importance. The rays with multiple reflections 
are considered to be weaker. We shall try to point out that 
this is not always valid and that also other than the rays 
with one reflection must be taken into account if we want to 
obtain a meaningful result. To demonstrate the nature of the 
approximation given by rays with only one reflection a pro- 
gram was written which adds all contributions of once- 
Bemiectéd rays (Fig. 6.5). 

The contribution pertinent to the ray reflected 


at the L-th interface is 


XG) = 18,, > %, 0 iwG(c) (6.34) 


Ls donk. fo) 1°12 n Met * 


where Xo is the amplitude at the point a . 


S = S = are transmission coefficients 


- Hoe be 2 oer ie Q,+Q,44 
across the i-th interface 
Q,-9 
Er ii =. Ltt is the coefficient of reflec- 
a eo 


tion from the L-th interface 
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Figure 6.5 


Plane waves propagating in the layered 
medium. 
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a = wits - an and K = pw = constant are vertical and 
horizontal wavenumbers in 
the i-th layer, respectively. 
G(ct) is a phase gain attributed to the ray path covered by 
the ray on its way down from point Gs to the reflecting 
interface and up from the reflecting interface to the 
receiver. 

The continuous velocity profile from Model IV was 
divided into homogeneous layers in the same manner as when 
the complete ray expansion was studied, i.e. into equal steps 
in velocity. The amplitude at the point Cs = t was set 


¥ 


equal to the exact value X, at - = Cy (20 km above the center 
of the transition). At each interface the contributions of 
rays reflected once at interfaces below were added to form 
the amplitude of the wave going upwards. Such addition 
depends on frequency and layer thicknesses and we observe an 
interference phenomenon. Thus the returning wave has varying 
amplitude which depends on frequency and the number of layers. 
Thus we cannot compare this value in any way to the complete 
ray expansion result or to the result for the continuous 
velocity profile. We can, however, compare the amplitudes at 
the turning point of the ray or of the downgoing wave which 


represents the first term in the ray series in Figure 6.4. The 


percentage error of this partial ray expansion 
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evaluated for the same parameters as in the previous example 
is displayed in Table 6.2 as a function of number of layers 
K . We see that the error does not decrease with the number 
of layers but remains rather large V14Z. *We%ean plot the 
absolute values of the percentage errors for the complete 
ray expansion and the partial ray expansion (Figure 6.6). 

The shaded area between the two curves can be interpreted as 
the difference due to the contribution of multiple reflected 
rays in this steady state case. The situation is slightly 
different if response to an impulsive source is investigated, 
especially for the first arrivals. Several authors showed 
that in this case once reflected rays give sufficiently 
accurate results (see, for example, Helmberger, 1968). For 
later arrivals the situation is different and multiply reflected 
rays must be taken into account. Muller (1970) computed 
theoretical seismograms for media consisting of homogeneous 
layers using rays with different number of reflections. He 
showed that only minor differences exist between the seismo- 
grams which include rays with only one reflection and those 
which include rays with three or five reflections, provided 


the medium has "moderate" number of layers (10-15 in his case). 
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THE PERCENTAGE ERROR OF THE PARTIAL RAY 
EXPANSION FOR MODEL IV AND w = 100 rad/sec 


Number of Layers Error at or 
K €[%] 
3 LS wey he 
4 - 3.618 
5 2.608 
7 - 4.987 
10 - 5.475 
LS = 10 se 
20 =10.15 
30 711.83 
40 mlz .69 
50 nh Mie Bea ap 
60 =-13.58 
70 $11./4 
80 -12.18 
90 -12.52 
100 -12.79 
120 =—13.20 
140 mA ae 
160 = 135 02 
180 =23593 
200 Bool Sete He Be 
250 algetz 
300 -13.49 
350 -13.76 
400 =i35..96 
450 -13.78 
500 -13.94 


Table 6.2 
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Nevertheless, for media with greater number of layers, the 
multiply reflected rays cannot be neglected. Waddington 
(1973) also showed that considerable differences exist 
between the seismograms which include only once reflected 
rays and those which include multiply reflected rays. 

The purpose of this study was to demonstrate the 
nature of the approximations of the solution for continuous 
velocity transition by complete or partial ray expansions. 

It was shown that the error in the complete ray expansion 
with respect to the exact solution is due to the bad 
approximation of the continuous model by layered medium when 
number of layers is small. The main source of the error in 
the partial ray expansion is in the choice of rays that take 
part in the partial ray expansion. With increasing number of 
layers the number of rays with three or five reflections 
grows considerably and their contribution plays a significant 


role in evaluation of the amplitude. 


6.3 Convergence of the ray expansion 


In any ray expansion the amplitudes of rays decrease 
with increasing number of reflections and refractions. This 
justifies the use of partial ray expansion in practise. The 
decrease in amplitude of one dynamic analogue may, however, 


be compensated by their number in the dynamic group. The 
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total dynamic effect of the group of dynamic analogues depends 


on 
of 
of 


by 


reflection and transmission coefficients 


are 


where 


In the case of 


for 


It is assumed that at the free 


near 


the dynamic analogue. 


the reflection coefficient 


their number in the group N 


In 


assuming only vertical or near vertical reflections. 


Ss NT 1 . 
re) 


q (6.22) and on the amplitude 


order to study the behaviour 


the total dynamic effect the problem will be simplified 


The 


at the i-th interface 


(5.39) 


reflections they become 





ll 
AG 


(6.36) 





surface z = Zo (Figure 6.1) 


29 The layers are 
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chosen in such a way that the reflection and transmission are 
constant (i.e. the indices in (6.36) can be omitted). Thus 
the amplitude of a dynamic analogue depends only on the total 
number of reflections and does not depend on their distribu- 
tion on the interfaces. The number of reflections at each 


interface is subject to conditions 


max(0,n,-n bi Se mes of n,-l Bett Lie dontaghi & C6..57 ) 


i+} i 


implied by (6.24). In analogy with the dynamic analogue 
characterized by number of reflections at each interface, 

the dynamic analogue characterized by total number of reflec- 
tions can be called a "generalized dynamic analogue". This 

is an important simplification which allows us to rewrite 

the total amplitude (6.24) of the group of kinematic analogues 


arriving at time (6.23) 


i=l ye ie 
in the following form 
M M 
max max 
E A, AY CNM) = & HL. (6.38) 
M=M m M=M ue 
min min 


/1¥ is the total number of reflections along the ray path 
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= : | 
= 2 ea | ears: Carey ot go ton) OM +620. = 1 (6.39) 
Es f=] i 1 a 


M is the total number of reflections from above the 


feeerteaces i @ l. «sa Jel 
M = >. mm (6.40) 


Every kinematic analogue (6.23) has total number of segments 


equal to 


J 
2N=2 Ef no (6.41) 


where n; is the half-number of segments in each layer. The 


minimum number of reflections from above the interfaces 


(6.42) 


M., = max(0,n,-n, 1) 


The maximum number of reflections from above the interfaces 


» * 1... seh d-lis,given.by 


J- 
We 3 Sims Be 
So J¢e] 


1 
(n,-1) =~ n=) = 2 (6.43) 


The total dynamic effect of groups of generalized dynamic 


analogues is denoted . my and the number of generaliZed ana- 







4 
f=L 
- i + | ’ ‘ v- a} as + »= 
L ou; ee) oe 
S) ods mo27? etolissltes to tedauvn Leios od? ef 
[- L ‘ee 7 f = t ‘bSoselts a 
L=T 
0 2 &€ # MM + 
> 4 
a iJ} : BI03 ‘2 ra | (é »* d) > vaolan & 52 Jemsata m 12) 
of 
. 
$ iA 
. 7 
b '. 
yf” se29 al siaemg: lo todmyna-tied sdi el 7 8% 
S J 
€ , 
vods mott enoizgoesfitiset ico tsdmun sueiegs 
: 7 
7 
vd qovig -ak I-l 2. . 
I=L , ‘ 
($d.8) (,,,0-,98,0)xem 2 = M 


“ irs 
7 = agankzerns ts 


: : 7 


svods mox? anolsasitez 20 tedaua auolxem 9) 
Far 7 a * ‘ - : mos 






o 
i] oo. 


7 re Te Peer 
iT Je w - d nev © pt [-t 


ey). i Pie’ 





183 


logues in the group is 7{(N,M). 


He. = baer’ ; 
tg AM ((N,M) (6.44) 
where 
Mtn _-1 Mm + 1e-n 
Rh), gavesijpie toiucidy J (6.45) 


2v 


Ii 
is the amplitude of the generalized dynamic analogue 
characterized by It reflections. Thus the total amplitude 
of all rays (or all kinematic analogues) arriving at time 

T is the sum (6.38) of total dynamic effects of groups of 
generalized dynamic analogues. 7 (Nn ,M) is the total number 
of generalized dynamic analogues in the group of generalized 


dynamic analogues with It cetlections: 


K 
M 
4 hi (N,M) = es N (9, >m, CK) (6.46) 


and the summation is over all possible distributions of M 
reflections into J-l layers so that in the i-th layer the 
number of reflections satisfies the condition (6.37) 

N (n, ,m, (k)) is the number of dynamic analogues in the group, 
where number of reflections m at i-th interface is given 

by k-th distribution of M . The number of distributions of 

M reflections into J-1 layers, subject to condition (6.37) is 


denoted by K To determine this we must solve the problem 


M . 


of determination of the number of distinct ways in which M 
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balls can be distributed into J-l different pockets allowing 


at most I, batten 1-th pocket (1 = 1, ..% J-1). 


5 ad n,-1 - Max(0,n, ) (6.47) 


i Pa ey 


Eisen (1969; p. 96) gives the solution to a simpler problem 
of distinct distributions of M balls into J-1l pockets allowing 
at most I balls per pocket. Then the number of these distri- 


butions is given by a simple recurrent formula 
r 
Ky = K(M; J-1,1)-=—2 — K(M-233-2,1) (6.48) 


and its results are displayed in Table 6.3 for the case of 
maximum number of balls allowed per pocket being I =3 . 


By convention 


ECO; O27). eat and E(IMsOe1) =.0 for M #0 
(6.49) 
The formula (6.48) can be generalized for the case, when the 


maximum number of balls allowed per pocket is different for 


different pockets. Generally, we can write 


pace te Vpeae (6.50) 


A 
I 


K(M;J-1,1) = K(M;I,,1 


rat Looe I Le a ene generate the recurrent 


formula analogously to that in (6.48): 
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THE NUMBER OF DISTINCT WAYS IN WHICH M BALLS CAN BE 
DISTRIBUTED INTO J-1 POCKETS ALLOWING AT 


MOST 3 BALLS PER POCKET 





Table 6.3 (from Eisen, 1969; p. 97) 
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fee. tt the first pocket contains £ balls, where 2 = 1 or 2 ...; 


or Il> the remaining M-& balls can be distributed into J-2 


“remaining pockets in K(M-231,, wk kins T5_y) Wave.  Tanhle 6.4 
gives the numbers of distributions for I, = 2, I, = 3, I, = 2, 
I, = 4 . The initial values remain the same as in (6.49). 

For every distribution of N, given by n,5i vg a ieee SA ae 


i.e. for every kinematic analogue, each of the distributions 
k, k=1, ... K, of M reflections represents a group of dynamic 


analogues N,(n,.m,(k)). Thus 


M 
x K,, = number of all groups of dynamic 
M=M analogues 


When investigating Table 6.4 we see, that every row has a 


maximum for 


M 
eee t (6.52) 


rounded up or down if M is odd. The fraction rounded down 


in denoted 


4 M 
Ms (= (6.53) 
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THE NUMBER OF DISTINCT WAYS IN WHICH M BALLS CAN BE 
DISTRIBUTED INTO J-1 POCKETS ALLOWING AT MOST 2 BALLS 
in TEE Biko. AND THIRD POCKETS, 3 BALLS IN THE 


SECOND POCKET AND 4 BALLS IN THE FOURTH POCKET 





Table) 64% 
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After studying the Table 6.3 it can be seen that the numbers 
of distribution K.., M = M ; 2 oe. 1M in each row must 
M min max 


always be smaller than the maximum in the previous row added 


E°+0l°times;: THOs (Sor kT care 


K We (ke, 6 (tian 7 


“ a ** ime Pe te (6.54) 


min max 


In the general case when the number of balls in the pockets 
is not equal for all pockets the numbers of distribution Ku 


can be bounded by 


[Max(I,,j=1, Dae Sen Giles Ss Ea err aaa 


A 
” 


(6.55) 


Every group of generalized dynamic analogues contains 
one group of dynamic analogues with maximum number of analogues 
for given M. We will denote this number of dynamic analogues 
by mg 687 omy (KE) and we can write an upper limit for TT oxy) 


as follows: 


2 


TE cum) << [Max(1,,J21, ... J-1)177? + ,(nysm,(k)) (6.56) 


j 


To find the "most powerful" group of dynamic analogues 
yg 695 °™y CK) is analytically very difficult since from (6.22) 


J-1 n n -1 
No(n,.mj(k)) = Tc i ¢ it3 
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J-1 

where M= fF m and m, (k) is given by the k-th distribution 
i=1 

of M, subject to the restriction (6.37). We were not able to 


find an analytic way of determining mg (Pym): Obviously 


an upper limit to (6.22) can be given by 


n n-L 
n oe ) (6.5) | 
2 


yN gq 694 om, Ck) < (Cc 
aed 


i} 
nae 
w 


h = M i -1). i i 
where n ax(n, ,i J-1) This overestimates mg 694 2™, 0K) 
but we do not know a better way at this moment. The number of 


generalized dynamic analogues in a group with Tk veflectious 


is bounded by 


1Y (NM) ¢ — CGnsb) G C ) (6.58) 


and total dynamic effect of the group of generalized dynamic 


analogues is 





Ane | < lay: —— ((n-1) c™ + c®t) (6.59) 


N 


Thus if we can prove that the upper limit of a tends to 


zero faster than (M —M_._) as n > ©, the convergence of 
max min 


the partial ray expansion has been proved. 


For large values of n the Stirling's approximation 


and ok fa TheS7). 


=-1 


(C.11) can be used to approximate co 
Zz 
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n+1 J-l1 

n n-l e ,2n J n+1 
“3 ae ~ {= 2 EL Sepa } for n even 

2 2 

(6.60) 

n n e) .2n 1 n ch See 
Cc Cc ~ {= 2 - = + (—~) } for n odd 

n n n n+1 

Caer hie 


~The only estimate we could find for the amplitude of one 


fxrgm €6.44) is 


generalized dynamic analog Aine 
por 1 n 2n_-n 
a gv JK 1 év bok BL 
beeen loge) <a (6.61) 


If we substitute (6.60) and (6.61) into (6.59), we see that 
the upper limit of Any | grows with n. Thus all the esti- 
Mates were too reugh to yield a useful result. 

Many attempts have been made to find a finer 
estimate and prove the convergence in the general case, but 
all with no success. The convergence was proved only for one 
special group of kinematic analogues which is the group of 
kinematic analogues with uniform distribution of segments into 
the layers (Fig. 6.7). The symmetry of the problem intro- 
duces symmetry into the "most powerful" group of dynamic 
analogues yg 693 9m KD) which, in the general case, had to 
be introduced artificially (see Eq. 6.57)... Thie, as we shall 
see, is an important property and allows us to give more 


realistic estimates to all the factors that contribute to the 
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Figure 6.7 Some of the dynamic analogues 
in the group of kinematic 
analogues Ni. (2,2,2). 
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total dynamic effect of the group of generalized dynamic 
analogues. 
We will denote n a half-number of segments in each 


layer. Thus 


nae on = 7, oor ai . i ae = (J-1) + (n-1) (6.62) 
and O< ms < n-l 
J-1 
a. M =~; = m ; then the recurrent formula (6.48) for the numbers 
i=1 


of distributions of M into J-1l layers is valid if I is re- 
placed by n-l. The "most powerful" group of dynamic analogues 
is that with number of reflections m being equal to oe 
(fraction rounded down). Thus the number of generalized 


analogues with hie reflections can be estimated by 


Pa a i 
TE-GUMY SS Kept oy Ny (5.05) = 8 (NM) (6.63) 
where 
+) Jui 
— e ,2n n-1'. 41." ": (6.64) 
Xn, 4) = ~——- ee 2 are Pier} } for n even 
n+] J~-1 
_— i * 1 e 2n n-l 
ea) ot es Pe Sree eet } for n odd 


After a considerable amount of algebra it can be shown that 
7¥ (N,M) has maximum at M and the function is diagramatically 


drawn in Figure 6.8. 
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Figure 6.8 





Schematic diagram of the behaviour of the ampli- 
tude [A my of generalized dynamic analogues, of 


the number of generalized dynamic analogues 
‘(N,M), and of the total dynamic effect : 
AM) An! 


The scales along the ordinate are different for 
all three functions. The only purpose is to show 
the approximate shape of those functions. 





194 


The amplitude of the generalized analogue charac- 


terized by 7? = ae + 2n - 1 reflections can be estimated for all 


n (both even and odd) as 


év 5 Sv, (n-1) (J-1) +1 


lAggel = Gy) aC (6.65) 


The upper estimate of the dynamic effect pertaining to the 


group of generalized dynamic analogues with ii reflections is 


a | < Ag > FRCN,M) (6.66) 


év 
and for ao < 45 we can write 


lim |fp-| = 0 (6.67) 
Rag ne 


fos] 


So far we have proved that the total dynamic effect pertinent 
to the most numerous or "most powerful" group of generalized 
dynamic analogues decreases with increasing half-number of 
segments n. But this group is only one of their total 

number Lae + 1 and represents only one term in the series 
(6.38). The next step is to show that any of the Te + ] 


groups of generalized dynamic analogues decreases with 


increasing n. If M is the point where | A has its maximum 


me 
a | ae | ated (6.68) 
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We are able to find its position with respect to M: 


-) d|A.,| | 
d 
ee Ts mw lige te rs C6469) 
where 
d|A,, | IN +1-n 
mu = SY) 0? log (32) 0 for M > 0 (6.70) 
and M <M 
SWS 0. for. Mo=-M (6.71) 
M>M 


re 


Therefore M < M (i.e. I” = 2(M4+n)-1 <M) and the behaviour of 
A | and 7¥(N,M) as functions of M is diagramatically drawn 
in Figure 6.8. The derivative 


dA yp Ris 


aM as Lo (6.72) 


and therefore M > M as n > © and the maximum of [AE ne | 
approaches to the maximum of J€(N,M). So we can say that 
for any € > 0 real, there exists My such that for any M > M, 


the following relation is valid 


IM ne | < lt aye | tné (6.73) 


where 7/{ = 2(M+n)-1 and IK = 2(M +n) -1 - Our objective is to 
investigate the convergence of the partial ray expansion when 


the amplitude of individual events is estimated by 
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max " 
U= re '(N,M) Nene J , # = 2M+ 2n- 1 
n m=o 8 67 
(6.74) 
If we use for the upper limit of an 
< ° ~ 
PaloSeMeae “dA He (6.75) 


where rs = (n-1)(J-1) and Tt = 2M + 2n - 1 we can prove that 


lim |u_| =0 
n 
n>o 
and therefore also 
lim. .U..+90 
n7~o n 


Thus we have proved the convergence of the partial ray series 
which consist only of kinematic analogues with uniform 
distribution of segments into individual layers. We hope 
that the procedure showed the enormous difficulties encoun- 
tered in the attempt to prove the convergence generally. It 
looks that the decision about significance of contribution 
from different groups of dynamic analogues must be done 
numerically for individual cases. We have not found any 
theoretical criterion according to which the ray series can 
be terminated at certain points with assurance that the 


remainder is negligable. 
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CONCLUSIONS 


In this thesis the non-geometrical effects of the 
wave propagation from a point source in the Epstein medium 
were studied. The response of the medium in the frequency 
domain was expressed as an integral in the complex plane of 
the horizontal wavenumber. The integrand contained the source 
strength and the reflection coefficient that describes the 
change of the wave on its path from the source to the receiver. 
The reflection coefficient for Epstein medium is an analytic 
function and the integration was performed both analytically 
and numerically. The analytic evaluation is an asymptotic 
one valid only at high frequencies. The numerical integration 
can be performed for any frequency needed for the determination 
of the spectrum used to synthesize the seismograms. The approach 
used here to study acoustic waves in fluid media is well 
suited to SH waves in elastic media. In principle it can also 
be applied to investigation of P-SV waves in elastic media 
provided a suitable potential representation is employed that 
results in the decoupling of P and SV waves at high frequencies. 

The behaviour of the response integral in the 
vicinity of caustics shows large amplitude decay at high 
frequencies. The amplitudes are much stronger at low frequencies 
which is the result of stronger partial reflection for longer 


wavelengths. In the region before the caustic that is the 
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illuminated region, the oscillation of amplitude due to the 
interference of two arriving waves is observed. The maximum 
amplitude is shifted from the geometrical endpoint into the 
illuminated region. For narrow transitions the behaviour near 
the caustic changes to critical point behaviour characterized 
by interference of the totally reflected wave and the head 
wave. The partial reflection is stronger and more like the 
sub-critical reflection from a sharp velocity increase. It is 
evident that the effect of finite dimensions of the transition 
layer on partially reflected wave and the head wave grows 

with increasing thickness. The dependence of the head wave 
amplitude changes with increasing thickness from = to 1 and 
should be useful for the interpretation of data when the 
transition thickness is determined. 

For transitions with velocity reversal the amplitude 
of the wave propagating along the reversal decays exponentially 
with distance. The “shadow" caused by this amplitude decrease 
is frequency independent which is quite different from a "true" 
shadow. The amplitude decay into the shadow caused by a 


lbh Therefore when 


sudden velocity decrease depends on exp (-aw 
interpretation is done these phenomena can be used to distin- 


guish between these two cases which on the travel-time curves 


have very similar characteristics. 
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The exact solutions for Epstein velocity profiles 
can be compared with approximate solutions for layered media. 
When the continuous Epstein structure is replaced by a layered 
medium, a procedure often used in seismology to find approxi- 
mate solutions for inhomogeneous media, the error arising 
from this approximation can be estimated. This approximation, 
which depends on the number of layers and the wavelength of 
propagating waves, is good only if a very large number of 
layers is used. With the computers available at the present 
time, it is impossible to compute the synthetic seismograms 
by Cagniard-de Hoop technique for that many layers. Therefore 
we must be satisfied with less accurate approximations of the 
medium. But even for a small number of layers further 
approximations are made by considering waves that contribute 
most to the total amplitude and neglecting the rest. The 
study of the partial ray expansion shows that not only the 
once reflected waves are significant but that multiply reflected 
waves often contribute significantly as well. Their contri- 
bution grows with growing number of layers and unless they 
are included in the evaluation of the amplitude the results 


are inconclusive. 
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APPENDIX A TRANSCENDENTAL FUNCTIONS 


This appendix contains details about some proper- 
ties of the hypergeometric and Airy functions frequently 
used in this thesis. The results are stated for reference 


purposes and their justification is not included. 


A.1l Hypergeometric functions 


The formulae in this section are taken from Erdelyi 
Seal. (1953, Chanter 2); Copson (1935, Chapter 10) and 
Abramowitz and Stegun CRS oa SO) s 


The hypergeometric function 





be (a) (db) ee 
F(a,b;c3&) = ¢£ Co) o) pen Awa! 
aot) € ns 
where 
ool (atn) 
SU es 
is a solution of the hypergeometric equation 
d*6 do 
§¢(14€) eo | {c-(a+b+1)—&} == - abo = 0 Le ee 
aE dé 


The equation A.1.2 has three regular singularities at 4 = Q, 
G5 = 1 and o3 = ©. Copson (1935, p. 237) shows that a second 
order differential equation with regular singularities has 


two independent solutions in terms of series. These series are 
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convergent in certain regions near the regular singularities. 


For equation A.1.2 the series solutions are: 


at €, =.0 o) = FF Cadopegs) 
l-c 
®, = (-&) F (a-c+1,b-c+1;2-c;&) 
ce ae ®, = F(a,b;atb+1l-c;1-é) 
2, = Gaye F (c-b ,c-a3c-a-bt+13;1-é) 
-b 1 
2c Se feo ®, = (-&) F(b,b-c+ls3b-at1;= ) 


uf 


ae F(a,a-ctl;a-b+1;& ) 


The hypergeometric function e = F(a,b;c;€) is defined by the 
series solution in region | = | < land by the analytic 
Esntinuation for |&| = 1. Tt ts regular for |&| < 1 and we 
shall see from its analytic continuation that it has two 
branch points at b5 = 1 and G3 =o , 

The analytic continuation is possible due to the 
fact that regions of convergence around each singularity have 
common regions through which the line of continuation can be 


drawn. At the same time, an integral solution was found by 
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valid everywhere in € plane with 


G5 and E43 


s 
(cts ie = 


The evaluation of this integral in different regions 


of € plane which correspond to the regions of convergence 


around the singularities gives the 


tion for every solution in A.1.3. 


correct analytic continua- 


By contour integration in 


the complex s plane along curve Cy (Fig. A.1) we get for 
be} <td 
aha Te? ae 
T, (&) = URE EX a. DsiGsc) Poca 
(Copson, 1935: p. 254). Cc, avoids all poles of the integrand 
in A.1.4 except those of I(-s). 


For |&| > 0 


is performed along curve C 


Thus we get 


T, (6) Pte-a) 


['(a-b) I (b) 


I'(c-b) 


_ I'(b-a)T (a) 


the integration in the complex s plane 


D which avoids all poles of I(-s). 
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Figure A.1 The complex s plane showing contour 
for Barnes integral. 
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for |E—| > 1. The formula 


as ['(a-b)T(c) x [(b-a)T(c) 6 
1 [(e-b)l (a) °5 P(c-a)I(b) 6 
gives the correct analytic continuation of 1: It defines 
the hypergeometric function o for |&| > 1. Similar formulae 
can be obtained for all hypergeometric functions A.1.3. 

The method of analytic continuation is well demon- 
strated through matrix notation. If 12 denotes a column 
vector of solutions convergent around o3 and iA a matrix 
of analytic continuation from region e. to region oy? we 


can write 


ie Yh? a.) = 1,2,3 ee ae 
o 6 o 
1 3 5 

Pie TE Fo ndpn td but( Mok gowgd ie Cd) Arals.9 
on ¢, — o, — o. 


The following expression demonstrates the analytic continuation 


of i2 throughout the entire € plane: 
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The elements of the 


It is easy to realize that 


13% 318° 


matrices of analytic continuation can be found from the 


connection formulae between different solutions which are 
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2, 


listed in Erdelyi et al. (1953, p. 106). We are interested 


mainly in matrix 314 which is 
A nN ['(1-c) I (1+b-a) '(c-1)T(1+b-a) 
a! 52 '(l-a)T(1+b-c) [(c-a)I(b) 
315 - 
A A (1l-c)P(1+a-b) ['(c-1)I (1+a-b) 





61 62 Pfiepyr(i+a=c) I (c-b)T (a) 


pA rad | 


The indices of the elements of the matrix are due to the 


components of the vectors they are connecting: 


6 = = = ,.A 9 sah a Be 


v6 461 = Ao oe) 


This may seem redundant but for use in Chapter 3 it proves 


useful. Similarly, for matrix 134 we have 





aera Aga “Aso 
1 
si 26 GA61 Asy ard 
where det A = a-b = (det A) 
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A.2 Airy functions and Stokes phenomenon 
oD LR ESSN &: NSN aera eelattetcRec Pel set aatcssheda cde asndlata 


The results and formulae presented in this appendix 
are taken from Budden (1961, Chapter 15) and Abramowitz and 
stegun (1965, p. 446). 


Solutions to the differential equation 


d*w 
dy? 


~ Sy woe. 0 > IPG 
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are Airy functions Ai(x), Bi(x) or Ai(x e ) or any linear 





combination. The following linear relation connects them 
2m 1 
Ai (x) Ai(xe 7) 
se = _2n4 Bade & 
Ba-(y) Ai(x e 3 ) 


The elements of the matrix are given by Abramowitz and Stegun 


K19655 pi 446) 


Ti mi 
3 = 
C= e 
M = Ter Fs 
7 ma _mi 
6 
e e 


Bot Teal oy <0. Ad (y} arid Bi(x) represent oscillatory solutions 
which in wave propagation can be thought of as standing waves. 
For xX > 0 real, Ai(x) is exponentially decaying while Bi(y) 


grows exponentially (Fig. A.2). This is evident from their 
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Figure A.2 





Airy functions Ai(x), Bi(x) for 
real values of xy . (Budden, 1961; 
pe 290). 
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asymptotic expressions which are in fact linear combinations 


of the WKBJ solutions to A.2.1: 





2 yol? 

SH) ab! 7 ay 
Ai(x) = kn? x es [arg x| <7 

-} - 2 rd T 2 
at Cy), = 1 i lx | % sin G lx | / + va larg (-x) | < 37 

Aw2.4 
+274 

From this we get asymptotic formulae for Ai(x e ’ ) on the 


negative real x axis in the form: 





Ai(x e a eat Et e are xX] < % 
AMZ. 5 
a 
Thus in application to wave propagation Ai(x e ) represents 
a wave travelling in positive y direction while Ai(x e 3 ) 


is a wave travelling in negative y direction. 

The asymptotic behaviour of the function Ai(x) is 
described in different regions of the X plane by different 
linear combinations of the WKBJ solutions. This property of 
the asymptotic approximations to the solutions of the differ- 
ential equations is called "Stokes phenomenon". The WKBJ 


solutions of A.2.1 are 
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They are in fact the second order saddle point contributions 


of the integrand in Airy integral tric. £.3) 


coo. 
1 Wea AS, 
Ai(x) = OTL e dy Ce a 
a 
Two saddle points at y = +x? give the contributions 

ae 2 

OI x" 

Let x e 


which differ from the WKBJ solutions A.2.6 only by constant 


an ~ = the exponents in 


A.2.6 are purely real. These lines are called the Stokes 


factors. On the lines arg yx = 0, 


lines and one solution always exponentially grows with lx | 
increasing and we call this the dominant solution. The other 
decays and is called the subdominant solution. On the lines 
arg X = 2. T, = x the exponent in A.2.6 is purely imaginary, 
the solutions have equal moduli and have an oscillatory 
character. We call these lines the anti-Stokes lines. The 
region between two anti-Stokes lines are called the Stokes 
regions where both solutions are present - the dominant and 
the subdominant. If we evaluate the Airy integral A.2.7 and 
larg x | < a the contour Cy Passes only through one saddle 


point y = x? (Fig. A.4.a). For [arg x| = PA the second saddle 
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point contribution must be considered too as Cc, passes through 
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Figure A.4 


yal}. 


a) 


b) 


Complex y plane with Stokes lines (S$) 

and anti-Stokes lines (A) of the Airy 
functions. The position of saddle 

points y., em) changes with arg y in 

tT Es gh Ee | ors show the path of 
integration (a) for arg ¥ = 0, th) for 
arg X = 7 (Budden, 1961; p. 303 and 1 
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both saddles y = + % (Fig. Al4.5)% ‘This is why the Stokes 
phenomenon occurs. On the Stokes lines the error introduced 
by the saddle point method of evaluation of the integral is 
greater than the subdominant term obtained from the contri- 
bution of the other saddle point. This is why we can neglect 
the subdominant term on the Stokes lines and express the 
asymptotic approximation only as one WKBJ solution. The 
detailed study of the Stokes phenomena for all four Airy 
functions can be found in Budden (1961, p. 294). The Stokes 
diagram in Fig. A.5 shows the asymptotic behaviour of Ai(y) 
and Bi(x) in dependence on arg X . For any argument a line 
inside the circle means the subdominant term is present and 


a line outside the circle means the dominant term is present. 
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(a) 





a ae | (b) 





Figure A.5 Stokes diagram for the functions 
(a) CAA (yd, Cb) Baty). (Budden, 1961; 
Pe 294). 





APPENDIX B VALIDITY OF THE REFLECTION COEFFICIENT 


The reflection coefficient for Epstein velocity 


transitions was found in Chapter 3. From (3.16) we get 


“e -(c-1) . 
ee} ['(c-1) 0 (1l-a)T (1+b-c) . FCa-ct1,b-c+l1;2-c;-e ~) y) 
4s R(e=ayT (GY —e pay 
) 


F(a,b3c3-e 


(B21) 


EOr 4 x peas <0 . This is only valid when the leading term 
in the hypergeometric series is considered i.e. far from the 
transition region. The remaining terms are dependent on 


frequency and thickness, since from (3.18) for © constant 


c-l = 2ioq, =-= i S cos Oo 
ye 
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Ae 


” 
atb-e = + (1 + 1607(q> - 4(q? + q2))) 


For monotonic velocity transition atb-c = +1 . First we show 
that the sign in the third equation is arbitrary and does not 


effect the results. We introduce new variables 


atb = 28 where B= et (B.2) 
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and the reflection coefficient becomes 


Pic=1) P(1-a~54%) TP (1-a-52%) ; 





RO? *) em *f eer tees 
VE P(1-c) P(-a+5+h) I (-a+5+4s) 


peeactpic weriga: 7. 
F(atsth , ~atsth3cs-e ah 


(haa) 


The term in the first curly brackets is evidently invariable 


to the change of sign. It remains to prove that 


ds 


F(g,h;2-c3;-e A E F(gt+tl1,h+1;2-c;-e ~) (B.4) 
F(-h,-g;c;-e ag F (-h+1,-g+l;c;-e a) 


From Abramowitz and Stegun (1965) formula 15.2.14 we get 


that the L.H.S. of the equation (B.4) is equal to 
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From formula 15.2.19 of Abramowitz and Stegun (1965) 
we get*that the RiH.S* of /(B.4)-is equni—to (€8.35)¢leThvs the 
expression (B.3) is invariable with respect to the charge of 
sign, To assess the validity of the reflection coefficient 
in (B.1) we must evaluate the absolute values of the second 


terms in the hypergeometric series: 
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F(a,bs;c;-e =) =1- _ Pwr ie o(e Gs, 
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for V4 < Vo * Thus the condition for validity of the reflec- 


> B.6 


Where 3 x "oe » the same condition applies to 4x concurrently. 


The transmission coefficient (3.17) for VB: = 4s is 
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The transmission coefficient is Ved 2 4 
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The condition (B.6) is satisfied for models (4.49), 
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this thesis, 
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APPENDIX C CONVERGENCE OF THE RESPONSE INTEGRALS 


To evaluate the response integral (2.22) we must 


evaluate the inverse Bessel transform 
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where R(T 5o 5K) Le piven by (4.2) or (5.2). If we use the 


following properties of the Bessel functions 
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TigC Kn) 3 [H, (kr) + E (xz).] 
aS?) (xr) - -H'7) (~cr) (c.2) 
the response integral (C.1) becomes 
+00 
P(w,r,0_) = - — P , (w) J R(G65o 5K) a? (er) a dk 
s s 
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The integrand in (C.3) is a multivalued function due to vertical 
wavenumber q . The complex kK plane consists of different 
Riemann surfaces given by different values of the square 

roots, The passage from one Riemann surface at a branch 


cut which is defined as any line connecting two branch points 
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often via infinity. Its choice is arbitrary, depending on 
the problem solved. The contour of integration along the 
real axis can be moved into any suitable position, provided 
it starts and ends on the same Riemann sheet, and, at the 
same time, the contributions from any singularities enclosed 
Or crossed are taken into account. 
For the velocity reversal the reflection coefficient 


per eaven by (5.2) 


P(1+2ioq,) P(s-2i0q,-y)T(4-2icgq +y ) wee Gn Somes on CY 
R(t a K) =< au net ey ac De 2 aaa aad *. -e ld See 
ia ofl P(1-2i0q,) lr (4-y) TP Cst+y) 


and the contour of integration in (C.3) can be moved to 
position in Figure 5.5 which is the steepest descent path. 
There are no poles in the first quadrant beyond the branch 

cut nor in the second. Only the poles + Sn (5.12). on the real 
axis exist and the contour is taken initially above these. 

The contribution along the arcs tends to zero as their 


radius grows to infinite values. For kK large the asymptotic 


expression for the reflection coefficient (5.22) may be used: 
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For large argument the asymptotic expression for aS) (xr) a 


used 


a tt 
(1) ikr-i Z 


1 2 
Be ete io (c.4) 





(formula 9.2.3 in Abramowitz and Stegun, 1965). 
(1) If kK = Rek + iImk and Rex > + © while ImkK = 0 we can 
do the following estimates: 
I——| v1 
qs 


aS”) (Kcr) ay Ik 7? + 0 (O25) 


and the reflection coefficient 


c, “im q + Im q 
HCE tk) | we © se Si eae (C.6) 


since Im qd, >"O,,' Tin qd = Or - < 0 and ce eS 


(2) If ImkK + + © and Rex > O 90r "Ret > s+ —| » the following 
1 
estimates can be made: 
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This proves that the contributions due to integration along 
the arcs can be neglected. Thus for the velocity reversal 
the main contribution to the response integral (C.3) is from 
the integration along the steepest descent path in Figures 
aoe ano 6.6. 

At high frequencies where the saddle point approxima- 
tion is valid the saddle point elias to — represents wave 
effected most by the velocity reversal are explained in 
Chapter 5. This saddle point approaches the pole Ky Ns Br Bs & 
(5.17) and is influenced by it. The integration contour can 
be moved so that it passes between the j-th and (ji+1)-th pole 
of the string ok, if we include their contribution. The 


residue contribution of the reflection coefficient at the 


j-th pole is 
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The residue series whose terms are (C.9) is not convergent. 
This becomes evident if the Stirling approximation for the 


gamma function is used 


1 — — 
T(z) % (27)? pe ‘6 foro are 2h at 


CG.14) 


(formula 6.1.37 in Abramowitz and Stegun, 1965). Then 
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because Im a4 aS) Tan qo Ce Ce mn Ce < 0... \Gnly.tf the 
exponential decay of the Hankel function due to distance over- 
comes the exponential growth of the reflection coefficient, the 
residue series will decrease for poles up to certain order. 

This happens at large distances r and then, only the zeroth pole 
is important and the rest can be neglected. It is not possible, 


however, to write the response as a complete residue series; 


we can only include the low order poles, say k, whose contribu- 
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tion decreases as order increases. The residue series forms 

an asymptotic series. That means the contour passes between 

k-th and (k+1)-th poles. The remainder of the contour intergration 
is of the same order of magnitude as the residue contribution 

of the (k+1l)-th pole. If the higher order poles were included 

the integral would increase in magnitude. The order k when 

the series should be terminated depends on distance: as range 
increases k also increases. This is due to the rapid decay 


given by distance between the poles 


r4.8 Im(,k KL) eS et Gee bes | 


k+1 2 


In such cases only the geroth order pole need be considered. 


If the residue of thezeroth order pole given by 
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(c.10), the residue of the first pole contribution 
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is greater and the residue series grows from the very beginning. 
This happens in the "illuminated" region and the integral 
cannot be evaluated by method of residues. 

For monotonic velocity transition the reflection 
coefficient in (C.3) is given by (4.2). The steepest descent 
path is shown in Figures 4.8 and 4.9. The integrand tends to 
zero in the second quadrant of complex kK plane and also in 
the first quadrant beyond the branch cut of qd, - This can 
be proved in the same way as that used above in the case of 
velocity peversed (Poe. (6.5), °(C.6), (C.7)} and (C.8)) provided 
we take the approximate reflection coefficient (4.17). The 
steepest descent path in Figure 4.8 ends in the valley above 
the real axis and if the original contour was deformed in 
this direction it would cross the branch cut. Thus we should 
return along the branch cut and then the contour can end 
anywhere in the fist quadrant beyond the branch cut. The 
steepest descent path in Figure 4.9 leads in the valley in 
the third quadrant. The contour cannot be distorted to end 
in the third quadrant as the integrand does not decay there 
with |k | growing. Again the integration should return to 


the real axis or first quadrant. 
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APPENDIX D ALGORITHMS FOR SPECIAL FUNCTIONS NEEDED FOR 
EVALUATION OF THE RESPONSE INTEGRAL 


D.1 Complex gamma function 


An algorithm for evaluation of the complex gamma 
function and its logarithm had to be developed in order to 
evaluate the reflection coefficient (3.22). At the time 
when this work had begun no standard program was available 
and we will describe here the algorithm developed for use 
in our computations, Since then, several algorithms appeared 
in literature as, for example, Algorithm 404, Complex gamma 
pumetton, by Lucas and Téerill (i971) and Algorithm 421, 
Complex gamma function with error control by Kuki (1972). 


The inverse gamma function and the logarithm 


8 
ae 2) 
log (T(z)) for z complex were computed using the following 


methods: 
CC®) +s DS €2) hi weg Ws 


(formula 6.1.23 from-Abramowitz and Stegun, 1965, henceforth 
referred to as Ref. 1) is used to restrict the argument into 


the first and second quadrants of the complex z plane 
l , 
T(z) = % T' Ze +11) (Ds daz) 


(6.1.15 in Ref. 1) is used to move the argument into the first 


quadrant so that 
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fa} RBe({z) «2.0 sf Inm(z2)' ¢ 2 


fb) TREC Z}esA10 ££ Im€z)or>«2 


The asymptotic formula is used to evaluate the log T(z) Et 


Re(z) 2 10 


1 1 
AE ba 1 cere 
log T(z) ~ (z-'5) log .z oe: slog (27) + oe ARE: 





(Duds) 


OG71F4Liin Referl)a olf tRe(z)r<+l0 othe following relations 
are used to restrict the argument so that 0 < Rete) © 1 and 


Beeeimiz) «< 0 . The Gauss multiplication formula 


log (nz) = 45(1-n) log, (27) + (nz-) log n + Log (2+) 


Ee eS 


Seat aciein Ref .°1) ts used if 1L.< Tn( 2) 6 2otosredice te 


0 < mor ye. 1 
T(z) = (2-1) (2-1) (D2 ooo 


(6.1.15 in Ref. 1) is used if Re(z) > 1; 0 < Im€z) < 1 to 
réduce to 0 < Re(z) < 1. Then Padé power series determined 


by Lee (1962) is used to evaluate the inverse gamma function 
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muvee nee) © 1, 0 < Im€z) << 1. This is too. accurate for 
an IBM 360/65 computer and some simpler series could be 
employed to improve the efficiency. 

The Fortran IV subroutines are used for real 
arguments after they were shifted to positive values 
performing step (D.1.2) above. 

When compared to Table 6.7 in Ref. 1 the results 


5 
are aocurate -to 1 part in 10 “ 


D.2 Hankel function of the zeroth order 


An algorithm for evaluation of the Hankel function 
aS) (2) for z complex was developed. It uses the Hankel's 
asymptotic expansion for | z | > > and the ascending series 
for }>| < 5. 

Rroem cermusa C9257), €922.9) and €9.2.10) in 
Abramowitz and Stegun (1965, from now on Ref. 1) we obtain 
the following asymptotic expansion for the Hankel function 


aS) (2) which is used if | z | 2s 
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Terms up to k = 3 are included in the computations. For small 
values of |z | the following method is used: 
Bere Ge te) ke dy fe) (D.2.4) 
re) Oo Oo a: 
where according to (9.1.12) and Coviels). in Ref. 2 
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yY is the Euler's constant. The first four terms are used for 
0 < [z| < 1.5 then next three terms are added LF 1°35 < jel < 3, 
and for 3 < |z| < 5 another three terms are added. 

When compared to Table 9.1 in Ref. 1 the results are 


accurate to 1 part in 10° : 
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